Method and apparatus for normalizing condition indicators

ABSTRACT

Disclosed are techniques used in connection with determining a health indicator (HI) of a component, such as that of an aircraft component. The HI is determined using condition indicators (CIs) which parameterize characteristics about a component minimizing possibility of a false alarm. Different algorithms are disclosed which may be used in determining one or more CIs. The HI may be determined using a normalized CI value. Techniques are also described in connection with selecting particular CIs that provide for maximizing separation between HI classifications. Given a particular HI at a point in time for a component, techniques are described for predicting a future state or health of the component using the Kalman filter. Techniques are described for estimating data values as an alternative to performing data acquisitions, as may be used when there is no pre-existing data.

CROSS REFERENCE TO RELATED APPLICATIONS

This application claims priority to U.S. Provisional patent applicationNo. 60/293,331 filed on May 24, 2001 which is incorporated by referenceherein.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This application relates to the field of vibration analysis and moreparticularly to performing vibration analysis for the purpose of devicemonitoring.

2. Description of Related Art

The transmission of power to rotors which propel helicopters and othershafts that propel devices within the aircraft induce vibrations in thesupporting structure. The vibrations occur at frequencies thatcorrespond to the shaft rotation rate, mesh rate, bearing passingfrequency, and harmonics thereof. The vibration is associated withtransmission error (TE). Increased levels of TE are associated withtransmission failure. Similar types of vibrations are produced bytransmissions in fixed installations as well.

Parts, such as those that may be included in a helicopter transmission,may be replaced in accordance with a predetermined maintenance and partsreplacement schedule. These schedules provide for replacement of partsprior to failure. The replacement schedules may indicate replacementtime intervals that are too aggressive resulting in needless replacementof working parts. This may result in incurring unnecessary costs asairplane parts are expensive. Additionally, new equipment may haveinstalled faulty or defective parts that may fail prematurely.

Thus it may be desirable to provide for an efficient technique fordetecting part and device degradation without unnecessarily replacingparts. It may be desirable that this technique also provide for problemdetermination and detection prior to failure.

SUMMARY OF THE INVENTION

In accordance with one aspect of the invention are a method executed ina computer system and computer program product for normalizing a set ofat least one observed condition indicator. A plurality of conditionalindicators and at least one associated factor are determined inaccordance with previous data acquisitions. A mean and at least onemodel coefficient corresponding to said at least one associated factorare determined. The set of observed condition indicators is adjusted inaccordance with model coefficients and said at least one associatedfactor producing a normalized set of condition indicators.

BRIEF DESCRIPTION OF DRAWINGS

Features and advantages of the present invention will become moreapparent from the following detailed description of exemplaryembodiments thereof taken in conjunction with the accompanying drawingsin which:

FIG. 1 is an example of an embodiment of a system that may be used inperforming vibration analysis and performing associated monitoringfunctions;

FIG. 2 is an example representation of a data structure that includesaircraft mechanical data;

FIG. 3 is an example of parameters that may be included in thetype-specific data portions when the descriptor type is an indexer;

FIG. 4 is an example of parameters that may be included in thetype-specific data portions when the descriptor type is anaccelerometer;

FIG. 5 is an example of parameters that may be included in thetype-specific data portions when the descriptor type is a shaft;

FIG. 6 is an example of parameters that may be included in thetype-specific data portions when the descriptor type is for a gear;

FIG. 7 is an example of parameters that may be included in thetype-specific data portions when the descriptor type is a planetarytype;

FIG. 8 is an example of parameters that may be included in thetype-specific data portions when the descriptor type is bearing type;

FIG. 9 is an example of a data structure that includes analysisinformation;

FIG. 10 is a more detailed example of an embodiment of a headerdescriptor of FIG. 9;

FIG. 11. is an example of a descriptor that may be included in theacquisition descriptor group of FIG. 9;

FIG. 12 is an example of a descriptor that may be included in theaccelerometer group of FIG. 9;

FIG. 13 is an example of a descriptor that may be included in the shaftdescriptor group of FIG. 9;

FIG. 14 is an example of a descriptor that may be included in the signalaverage descriptor group of FIG. 9;

FIG. 15 is an example of a descriptor that may be included in theenvelope descriptor group of FIG. 9;

FIG. 16 is an example of a planetary gear arrangement;

FIG. 17A is an example of an embodiment of a bearing;

FIG. 17B is an example of a cut along a line of FIG. 17A;

FIG. 18A is an example of a representation of data flow in vectortransformations;

FIG. 18B is an example of a representation of some of the CI algorithmsthat may be included in an embodiment, and some of the various inputsand outputs of each;

FIG. 19 is an example of a graphical representation of a probabilitydistribution function (PDF) of observed data;

FIG. 20 is an example of a graphical representation of a cumulativedistribution function (CDF) observed data following a gamma (5,20)distribution and the normal CDF;

FIG. 21 is an example of a graphical representation of the differencebetween the two CDFs of FIG. 20;

FIG. 22 is an example of a graphical representation of the PDF ofobserved data following a Gamma (5,20) distribution and a PDF of thenormal distribution;

FIG. 23 is an example of another graphical representation of the twoPDFs from FIG. 22 shown which quantities as intervals rather thancontinuous lines;

FIG. 24A is an example of a graphical representation of the differencesbetween the two PDFs of observed data and the normally distributed PDF;

FIGS. 24B-24D are examples of a graphical data displays in connectionwith a healthy system;

FIGS. 24E-24G are examples of graphical data displays in connection witha system having a fault;

FIG. 25 is a flowchart of steps of one embodiment for determining healthindicators (HIs);

FIG. 26 is a graphical illustration of the probability of a false alarm(PFA) in one example;

FIG. 27 is a graphical illustration of the probability of detection (PD)in one example;

FIG. 28 is a graphical illustration of the relationship between PD andPFA and threshold values in one embodiment;

FIG. 29 is an graphical illustration of the probability of Ho andthreshold values in one embodiment;

FIG. 30 is an example of an embodiment of a gear model;

FIG. 31 is a graphical representation of an estimated signal having aninner bearing fault; and

FIG. 32 is a graphical representation of the signal of FIG. 31 as afrequency spectrum.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT(S)

Referring now to FIG. 1, shown is an example of an embodiment of asystem 10 that may be used in performing vibration analysis andmonitoring of a machine such as a portion of an aircraft. The machinebeing monitored 12 may be a particular element within an aircraft.Sensors 14 a through 14 c are located on the machine to gather data fromone or more components of the machine. Data may be collected by thesensors 14 a through 14 c and sent to a processor or a VPU16 for datagathering and analysis. The VPU16 analyzes and gathers the data from theSensors 14 a through 14 c.

The VPU16 may also use other data in performing analysis. For example,the VPU16 may use collected data 18. One or more of the Algorithms 20may be used as input into the VPU16 in connection with analyzing datasuch as may be gathered from the Sensors 14 a through 14 c.Additionally, configuration data 22 may be used by the VPU16 inconnection with performing an analysis of the data received for examplefrom the Sensors 14 a through 14 c. Generally, configuration data mayinclude parameters and the like that may be stored in a configurationdata file. Each of these will be described in more detail in paragraphsthat follow.

The VPU16 may use as input the collected data 18, one or more of thealgorithms 20, and configuration data 22 to determine one or morecondition indicators or CIs. In turn, these condition indicators may beused in determining health indicators or HIs that may be stored forexample in CI and HI storage 28. CIs describe aspects about a particularcomponent that may be useful in making a determination about the stateor health of a component as may be reflected in an HI depending on oneor more CIs. Generally, as will be described in more detail inparagraphs that follow, CIs and HIs may be used in connection withdifferent techniques in determining an indication about monitoredcomponents such as Machine 12. As described in more detail elsewhereherein, the configuration data may include values for parameters thatmay vary in accordance with the type of the component being monitored.

It should be noted that the collected data 18 may include data collectedover a period of time from sensors such as 14 a through 14 c mounted onMachine 12. A user, such as a Pilot 26, may use a special serviceprocessor, such as the PPU24, connected to the Machine 12 to obtaindifferent types of data such as the CI and HI values 28.

As described in connection with FIG. 1, the VPU16 may receive inputsfrom Sensors 14 a through 14 c. These sensors may be different types ofdata gathering monitoring equipment including, for example, highresolution accelerometers and index sensors (indexors) or tachometersthat may be mounted on a component of Machine 12 at carefully selectedlocations throughout an aircraft. Data from these sensors may be sampledat high rates, for example, up to 100 kilohertz, in order for the VPU16to produce the necessary CI and HI indicators. Data from these sensorsand accelerometers may be acquired synchronously at precise intervals inmeasuring vibration and rotational speeds.

Generally, the different types of data gathering equipment such as 14a-14 c may be sensors or tachometers and accelerometers. Accelerometersmay provide instantaneous acceleration data along whatever axis on whichthey are mounted of a particular device.20 Accelerometers may be used ingathering vibration analysis data and accordingly may be positioned tooptimally monitor vibration generated by one or more mechanicalcomponents such as gears, shafts, bearings or planetary systems. Eachcomponent being monitored may generally be monitored using twoindependent sensors to provide confirmation of component faults and toenable detection of sensor faults.

No accelerometer is completely isolated from any other component. Thus,the component rotational frequencies share as few common divisors aspossible in order to maximize the effectiveness of the monitoringfunction being performed. For example, all gears being monitored shouldhave differing number of teeth and all bearings should have differingnumbers and sizes of balls or rollers. This may allow individualcomponents to be spectrally isolated from each other to the extent thattheir rotational frequencies are unique.

The indexers (index sensors) or tachometers may also be used as aparticular monitoring component 14 a through 14 c to gather data about aparticular component of Machine 12. The indexers produce a periodicanalog signal whose frequency is an integer multiple of theinstantaneous rotation frequency of the shaft that they are monitoring.These signals may be generated magnetically using one or more evenlyspaced metallic protrusions on the shaft passing by the fixed sensor.Alternatively, these may be monitored optically using a piece ofoptically reflective material affixed to the shaft. It should be notedthat each index point should be fixed in time as precisely as possible.In connection with magnetic sensors, this may be accomplished forexample by interpolating the zero crossing times of each index pulse andsimilarly for optical sensors by locating either rising or fallingedges. Assuming the minimal play or strain in the drive train whensomething is under load, the relative position and rate of any componentmay be calculated using a single index or wave form.

Because of the high data rates and lengthy processing intervals,diagnostics may be performed, for example, on pilot command or on apredetermined flight regime or time interval.

Each of the algorithms 20 produces one or more CIs described elsewhereherein in more detail. Generally, the CI may yield useful informationabout the health of a monitored component. This condition indicator orCI as well as HI may be used in determining or predicting faults ofdifferent components.

It should be noted that the VPU16 is intended to be used in a widevariety of mechanical and electrical environments. As described herein,different components of an aircraft may be monitored. However, this isonly one example of a type of environment in which the system describedherein may be used. As known to those skilled in the art, the generalprinciples and techniques described herein have much broader and generalapplicability beyond a specific aircraft environment that may used in anexample here.

In connection with the use of CIs, the VPU16 uses the CIs as input andportions of the data such as, for example, used in connection with analgorithm to provide HIs. These are described in more detail inparagraphs that follow.

It should be noted that in a particular embodiment, each mechanical partbeing monitored may have one or more sensors associated with it where asensor may include for example an accelerometer or a tachometer.Generally, accelerometers may be used, for example, to obtain dataregarding vibrations and a tachometer may be used, for example, to gaininformation and data regarding rotation or speed of a particular object.Data may be obtained and converted from the time to the frequencydomain.

A particular algorithm may provide one or more CIs. Each of thealgorithms may produce or be associated with a particular CI. One ormore CIs may be used in combination with a function to produce an HI fora particular part or type. As will be described in more detail herein,each of the algorithms may be associated or classified with a particularpart or type. The CI generally measures vibrations and applies afunction as described in accordance for each algorithm. Generally,vibration is a function of the rotational frequency in the amount oftorque. Using torque and a particular frequency, a CI is appropriatelydetermined in accordance with a selected algorithm for a part.

The algorithms 20 may be classified into four families or groups inaccordance with the different types of parts. In this example, thefamilies of algorithms may include shaft, gears, bearings, and planetarygears. Associated with each particular part being monitored may be anumber of CIs. Each CI may be the result or output of applying adifferent one of the algorithms for a particular family. For example, inone embodiment, each gear may have an associated 27 CIs, each bearingmay have 19 CIs, each shaft may have 22 CIs, and each planetary gear mayhave two or three CIs. It should be noted that each one of these numbersrepresents in this example a maximum number of CIs that may be used orassociated with a particular type in accordance with the number ofalgorithms associated with a particular class or family. Generally, thedifferent number of CIs that may be associated with a particular typesuch as a gear try to take into account the many different ways in whicha particular gear may fail. Thus, a CI reflects a particular aspect orcharacteristic about a gear with regard to how it may fail. Differenttechniques used in computing CIs are described, for example, in“Introduction to Machinery Analysis and Monitoring, Second Edition”,1993, Penn Well Publishing Company of Tulsa, Okla., ISBN 0-87814-401-3,and “Machinery Vibration: measurement and analysis”, 1991, McGraw-HillPublishing, ISBN-0-07-071936-5.

Referring now to FIG. 2, shown is an example of a data structure 50 thatincludes aircraft mechanical data. Generally, this data structureincludes one or more descriptors 56 a through 56 n. In this embodimentthere may be one descriptor for each sensor. A descriptor associatedwith a particular sensor includes the parameters relevant to theparticular component being monitored. Each of the descriptors such as 56a includes three portions of data. The field 52 identifies a particulartype of descriptor. Each of the descriptors also includes a common dataportion 54 which includes those data fields common to all descriptortypes. Also included is a type specific data portion 56 which includesdifferent data fields, for example, that may vary in accordance with thedescriptor type 52.

Descriptor types may include, for example, an indexer, an accelerometer,a shaft, a gear, a planetary gear, or a bearing descriptor type valuecorresponding to each of the different types of descriptors. The commondata portion 54 may include, for example, a name, part number andidentifier. In this example, the identifier in the common data filed 54may uniquely identify the component and type.

Referring now to FIGS. 3 through 8, what will be described are examplesof descriptor type specific parameters or information that may beincluded in a descriptor of a particular type, such as in area 56 of thedata structure 50.

Referring now to FIG. 3, shown is an example of parameters that may beincluded in a descriptor 60 which is an indexer descriptor type. Theparameters that may be included are a channel 62, a type 64, a shaftidentifier 66, a pulses per revolution parameter 68, a pulse widthparameter 70, and a frequency of interest 72 for this particular type ofdescriptor. It should be noted that the type in this example for theindex or descriptor may be one of sinusoidal, pulse such as 1/rev, oroptical. The shaft identifier 66 is that as may be read or viewed by theindexer that calculates the shaft rate. The pulse width 70 is in secondsas the unit value. Additionally, the frequency of interest 72 for thisdescriptor type is a nominal pulse frequency that is used in computingthe data quality signal to noise ratio. The use of these particular datastructures and parameters is described in more detail in paragraphs thatfollow.

Referring now to FIG. 4, shown is an example of the parameters that maybe included in an accelerometer descriptor type 80. The descriptor foran accelerometer type may include the channel 82, a type 84, asensitivity 86 and a frequency of interest 88. In this example for theaccelerometer descriptor type, the type may be one of normal, or remotecharge coupled. The frequency of interest may be used in computing thedata quality signal to noise ratio. The frequency of interest for a gearis the mesh rate which may be calculated from the gear shaft rate andthe number of teeth of the gear.

Referring now to FIG. 5, shown is an example of descriptor type specificparameters or data that may be included when a descriptor type is theshaft descriptor. A shaft descriptor 90 includes path parameter or data92 and nominal RPM data 94. The path data is an even length sequence ofgear tooth counts in the mechanical path between the shaft in questionand a reference shaft. The driving gears alternate with driven gearssuch that the expected frequency of a gear, shaft, bearing and the likemay be determined based on an input shaft RPM.

Referring now to FIG. 6, shown is an example of data or parameters thatmay be included in a descriptor when the descriptor type is the geardescriptor. Included in the gear descriptor 100 is the shaft identifier102 to which the gear is mounted and a parameter 104 indicating thenumber of teeth in the gear.

Referring now to FIG. 7, shown is an example of an embodiment of aplanetary descriptor 110 identifying those parameters or data that maybe included when the type is a planetary descriptor type. The planetarydescriptor 110 may include an input shaft identifier 112, an outputshaft identifier 114, a parameter indicating the number of planet gears116, a parameter indicating the number of teeth on the planet gear, aparameter 120 indicating the number of teeth on the ring gear, and aparameter 122 indicating the number of teeth on the sun gear. It shouldbe noted that the number of teeth on a planet gear relates to a planetcarrier that is assumed to be mounted to the output shaft. Additionally,the ring gear is described by parameter 120 is assumed to be stationeryand the sun gear 122 as related to parameter 122 is assumed to bemounted to the input shaft. It should be noted that the path between theinput and the output shaft may be reduced to using a value S for thedriving path tooth count and R+S as the driven path tooth count where Rand S are the ring and sun tooth counts respectively. An example of aplanetary type gear is described in more detail elsewhere herein.

Referring now to FIG. 8, shown is an example of a bearing descriptor130. The bearing descriptor 130 may include descriptor type specificfields including a shaft identifier 132, a cage ratio 134, a ball spinratio 136, an outer race ratio 138 and an inner race ratio 140. Anexample of a bearing is described in more detail elsewhere herein.

It should be noted that the data structures described in connection withFIGS. 2 through 8 are those that may be used in storing data obtainedand gathered by a sensor such as 14 a when monitoring a particularcomponent of a machine 12. Data may be gathered and stored in the datastructure for a particular descriptor or descriptors and sent to the VPU16 for processing. It should be noted that a particular set of data maybe gathered at a particular instance and time, for example, inconnection with the synchronous data gathering described elsewhereherein. In connection with this, a data set may include multipledescriptors from sampling data at a particular point in time which issent to the VPU 16.

What will now be described are those data structures that may beassociated with an analysis definition that consists of a specific dataacquisition and a subsequent processing of this data to produce a set ofindicators for each of the desired components.

Referring now to FIG. 9, shown is an example of the data structure 150that contains analysis data. Each instance of analysis data 150 asrepresented in the data structure includes a header descriptor 152 anddescriptor groups noted as 164. In this example there are fivedescriptor groups although the particular number may vary in anembodiment. Each of the descriptor groups 154 through 162 as identifiedby the group identifier 164 includes one or more descriptors associatedwith a particular group type. For example, descriptor group 154 is theacquisition group that includes a descriptor for each sensor to beacquired. The accelerometer group 156 consists of a descriptor for eachaccelerometer to be processed. The shaft group 158 includes a descriptorfor each shaft to be processed. The signal average group 160 includes adescriptor for each unique parameter set. The envelope group 162includes a descriptor for each unique parameter.

Referring now to FIG. 10, shown is a more detailed example of a headerdescriptor 170. Parameters that may be included in a header descriptor170 include: an analysis identifier 172, acquisition time out parameter174 and processing time out parameter 176. In this example, theacquisition, time out and processing time out parameters are in seconds.

Referring now to FIG. 11, shown is an example of a descriptor that maybe included in the acquisition group. A descriptor 180 included in theacquisition group may include a sensor identifier 182, a sample rateparameter in Hz 184, a sample duration in seconds 186, a gain controlsetting, such as “auto” or “fixed” 188, an automatic gain control (AGC)acquisition time in seconds 190, an automatic gain control (AGC)headroom factor as a number of bits 192 and a DC offset compensationenable 194.

Referring now to FIG. 12, shown is an example of a descriptor 200 thatmay be included in the accelerometer group. A descriptor in theaccelerometer group may include a parameter that is an accelerometeracquisition analysis group identifier 202, a list of associatedplanetary identifiers to be processed 204, a list of associated shaftanalysis group identifiers to be processed 206, a processor identifier208, a transient detection block size 210, a transient detection RMSfactor 212, a power spectrum decimation factor 214 specified as a powerof 2 and a power spectrum block size also specified as a power of 2.

In one embodiment, the list of associated planetary identifiers 204 alsoincludes two signal average analysis group identifiers for eachplanetary identifier, first identifier corresponding to the input shaftand a second corresponding to an output shaft.

It should be noted that the processor identifier 208 will be used inconnection with assigning processing to a particular DSP or digitalsignal processor.

Referring now to FIG. 13, shown is an example of an embodiment of adescriptor 280 that may be included in the shaft group. The descriptor220 may include a shaft identifier 222, a signal average analysis groupidentifier 224, a list of gear identifiers to be processed 226, a listof bearing identifiers to be processed 228 and a list of associatedenvelope analysis group identifiers 230.

Referring now to FIG. 14, shown is an example of a descriptor 232 thatmay be included in the signal average group. It should be noted that thesignal average group includes a descriptor for each unique parameterset. The signal average processing group is run for each accelerometerand shaft combination even if it has the same parameters as anothercombination. Each descriptor 232 may include a number of output pointsper revolution 234 and a number of revolutions to average 236.

Referring now to FIG. 15, shown is an example of a descriptor 240 thatmay included in the envelope group. It should be noted that the envelopegroup includes a descriptor for each unique parameter. It is notnecessary to repeat an envelope processing for each bearing if theparameters are the same. Each descriptor 240 may include a durationparameter 242 specifying the seconds of raw data to process, an FFT size244 which is a power of 2, a lower bound frequency in Hz 246, and anupper bound frequency, also in, Hz 248.

Referring now to FIG. 16, shown is an example of an embodiment 300 of aplanetary gear arrangement. Generally, a planetary gear arrangement asdescribed in connection with the different types of gears and items tobe monitored by the system 10 of FIG. 1 may include a plurality of gearsas configured, for example, in the embodiment 300. Included in thearrangement 300 is a ring gear 302 a plurality of planet gears 304 athrough 304 c and of sun gear 306. Generally, the gears that aredesignated as planets move around the sun gear similar to that as asolar system, hence the name of planet gear versus sun gear. Thearrangement shown in FIG. 16 is a downward view representing thedifferent types of gears included in an arrangement 300.

Referring now to FIG. 17A, shown is an example of an embodiment 320 of abearing. The bearing 320 includes a ring or track having one or morespherical or cylindrical elements (rolling elements) 324 moving in thedirection of circular rotation as indicated by the arrows. Differentcharacteristics about such a structure of a bearing may be important asdescribed in connection with this embodiment. One characteristic is an“inner race” which represents the circumference of circle 322 a of theinner portion of the ring. Similarly, the “outer race” or circumference322 b representing the outer portion of the ring may be a considerationin connection with a bearing.

Referring now to FIG. 17B, shown is an example of a cut along line 17Bof FIG. 17A. Generally, this is cut through the ring or track withinwhich a bearing or bearings 324 rotate in a circular direction. The ballbearings move in unison with respect to the shaft within a cage thatfollows a track as well as rotate around each of their own axis.

Referring now to FIG. 18A, shown is an example of a representation 550of different transformations that may be performed and the associateddata flow and dependencies for each particular sensor. The output of thetransformations are transformation vectors and may be used in additionto analysis data or raw data, such as bearing frequency, mesh frequency,and the like, by an algorithm in producing a CI.

Referring to the representation 550, an in going arrow represents dataflow input to a transformation. For example, the FF or Fast Fouriertransform takes as an input data from the A1 signal average datatransform. A1 has as input the accelerometer data AD. It should be notedthat other embodiments may produce different vectors and organize datainputs/outputs and intermediate calculations in a variety of differentways as known to those skilled in the art.

Referring now to FIG. 18B, shown is an example of a representation 350relating algorithms, a portion of input data, such as sometransformation vectors, and CIs produced for each type of component,that may be included in an embodiment. Other embodiments may usedifferent data entities in addition to those shown in connection withFIG. 18B. As described elsewhere herein, each type of component in thisexample is one of: indexer, accelerometer, shaft, gear, planetary, orbearing. Certain algorithms may be used in connection with determiningone or more CIs for more than one component type. It should be notedthat a variety of different algorithms may be used and are known by oneof ordinary skill in the art, as described elsewhere herein in moredetail. The following are examples of some of the different techniquesthat may be used in producing CIs. Additionally, FIG. 18B illustrates anexample of relationships between some algorithms, a portion of theirrespective inputs and outputs, as well as how the algorithms may beassociated with different component types. However, it should be notedthat this illustration is not all inclusive of all algorithms, allrespective inputs and outputs, and all component types.

What will now be described are algorithms and the one or more CIsproduced that may be included in an embodiment. It should be noted thatthe number and type of algorithms included may vary in accordance withan embodiment. Additionally, it should be noted that FIG. 18B may notinclude each and every input and output for an algorithm as describedherein and other embodiments of the algorithms described generallyherein may also vary.

The data quality (DQ) algorithm 356 may be used as a quality assurancetool for the DTD CI. DQ performs an assessment of the raw uncalibratedsensor data to insure that the entire system is performing nominally. DQmay be used to identify, for example, bad wiring connections, faultysensors, clipping, and other typical data acquisition problems. The DQindicator checks the output of an accelerometer for “bad data”. Such“bad data” causes the SI to be also be “bad” and should not be used indetermining health calculations.

What will now be described are the different indicators that may beincluded in an embodiment of the DQ algorithm. ADC Bit Use measures thenumber of ADC bits used in the current acquisition. The ADC board istypically a 16 bit processor. The log base 2 value of the maximum rawdata bit acquired is rounded up to the next highest integer. Channelswith inadequate dynamic range typically use less than 6 bits torepresent the entire dynamic range. ADC Sensor Range is the maximumrange of the raw acquired data. This range cannot exceed the operationalrange of the ADC board, and the threshold value of 32500 is just belowthe maximum permissible value of +32767 or −32768 when the absolutevalue is taken. Dynamic Range is similar to the ADC Sensor Range, exceptthe indicator reports dynamic channel range as a percent rather than afixed bit number. Clipping indicates the number of observations ofclipping in the raw data. For a specific gain value, the raw ADC bitvalues cannot exceed a specific calculated value. Low Frequency Slope(LowFreqSlope) and Low Frequency Intercept (lowFreqInt) use the first 10points of the power spectral density calculated from the raw data andperform a simple linear regression to obtain the intercept and slope inthe frequency-amplitude domain. SNR is the signal to noise ratioobserved in each specific data channel. A power spectral density iscalculated from the raw uncalibrated vibration data. For each datachannel, there are known frequencies associated with certain components.Examples include, but are not limited to, gear mesh frequencies, shaftrotation rates, and indexer pulse rates. SNR measures the rise of aknown tone (corrected for operational speed differences) above thetypical minimum baseline levels in a user-defined bandwidth (generally+/−8 bins).

The Statistics (ST) algorithm 360 is associated with producing aplurality of statistical indicators 360 a. The Root-Mean-Square (RMS)value of the raw vibration amplitude represents the overall energy levelof the vibration. The RMS value can be used to detect major overallchanges in the vibration level. The Peak-To-Peak value of the rawvibrating amplitude represents the difference between the two vibrationextrema. When failures occur, the vibration amplitude tends to increasein both upward and downward directions and thus the Peak-To-Peak valueincreases. The Skewness coefficient (which is the third statisticalmoment) measures the asymmetry of the probability density function(p.d.f.) of the raw vibration amplitude. Since it is generally believedthat the p.d.f. is near Gaussian and has a Skewness coefficient of zero,any large deviations of this value from zero may be an indication offaults. A localized defect in a machine usually results in impulsivepeaks in the raw vibration signal, which affects the tails of the p.d.f.of the vibration amplitude. The fourth moment (Kurtosis) of thedistribution has the ability to enhance the sensitivity of such tailchanges. It has a value of 3 (Gaussian distribution) when the machineryis healthy. Kurtosis values larger than 3.5 are usually an indication oflocalized defects. However, distributed defects such as wear tend tosmooth the distribution and thus decrease the Kurtosis values.

The ST algorithm may be performed on the following vectors: AD rawaccelerometer data, A1signal average data, RS residual data, NB narrowband data, and EV envelope data and others, some of which are listed in360 b.

The Tone andBase Energy algorithm(TB) 362 uses tone energy and baseenergy. Tone Energy is calculated as the sum of all the strong tones inthe raw vibration spectrum. Localized defects tend to increase theenergy levels of the strong tones. This indicator is designed to providean overall indication of localized defects. “Strong tones” aredetermined by applying a threshold which is set based on the mean of allthe energy contents in the spectrum. Any tones that are above thisthreshold are attributed to this indicator. The Base Energy measures theremaining energy level when all the strong tones are removed from theraw vibration spectrum. Certain failures such as wear, do not seem toaffect the strong tones created by shaft rotation and gear mesh, theenergy in the base of the spectrum could potentially be a powerfuldetection indicator for wear-related failures. Note that the sum of ToneEnergy and Base Energy equals the overall energy level in the spectrum.

SI are miscellaneous shaft indicators. SO1 (Shaft Order 1 in g) is theonce-per-rev energy in the signal average, and is used to detect shaftimbalance. SO2 (Shaft Order 2 in g) is the twice-per-rev energy in thesignal average, and is used to detect shaft misalignment. GDF (Geardetector fault) may be an effective detector for distributed gear faultssuch as wear and multiple tooth cracks, and is a complement of theindicator signalAverageL1 (also known as gearLocalFault).

In addition to the specifically referenced vectors below, the SIalgorithm takes input from the indexer zero-crossing vector (ZC).

The Demodulation analysis (DM) 370 is designed to further reveal sideband modulation by using the Hilbert transform on either the narrow bandsignal (narrow band demodulation) or the signal average itself (wideband demodulation) to produce the Amplitude Modulation (AM) and PhaseModulation (FM) signals. The procedures involved to obtain such signalsare:

Perform Hilbert transform on the narrow band signal (or signal average).

Compute the amplitude of the obtained complex analytic signal to obtainthe AM signal.

Compute the phase angles of the analytic signal to obtain the FM signal.

Compute the instantaneous amplitude of the analytic signal to obtain thedAM signal.

Compute the instantaneous phase angles of the analytic signal to obtainthe dFM signal.

The DM algorithm is performed on the band passed filtered data at afrequency of interest by taking a Hilbert Window function of thefrequency domain data and converting the data back to the time domain.

The Sideband Modulation (SM) 368 analysis is designed to reveal anysideband activities that may be the results of certain gear faults suchas eccentricity, misalignment, or looseness. CIs included in 368 a areDSMn. DSMn is an indicator that characterizes the Degree of SidebandModulation for the nth sideband (n=1, 2, and 3). The DSMn is calculatedas the sum of both the nth high and low sideband energies around thestrongest gear meshing harmonic. As indicated in 368 b, the SM algorithmis performed on the Fast Fourier transform vector (FF).

The Planetary Analysis (PL) 364 extracts the Amplitude Modulation (AM)signal produced by individual planet gears and compares the “uniformity”of all the modulation signals. In general, when each planet gear orbitsbetween the sun and the ring gears, its vibration modulates thevibration generated by the two gears. It is believed that when one ofthe planet gears is faulty, the amplitude modulation of that planet gearwould behave differently than the rest of the planet gears. Theprocedure to perform this algorithm is to obtain signal averages for theinput, output, and planet shafts. For each signal average:

Locate the strongest gear meshing harmonic.

Bandpass filter the signal average around this frequency, with thebandwidth equals to twice the number of planet gears.

Hilbert transform the bandpass filtered signal to obtain the AM signal.

Find the maximum(MAX) and minimum(MIN) of the AM signal.

Calculate the Planet Gear Fault (PGF) indicator as included in 364 aaccording to the equation PGF=MAX(AM)/MIN(AM).

The inputs to the PL algorithm are the raw accelerometer data (AD) andthe indexer zero-crossing data (ZC).

The Zero-Crossing Indicators (ZI) algorithm 354 is performed on thezero-crossing vector (ZC). The zero crossing indicators may bedetermined as follows:

D_(j)=In_(j+1)−In_(j), j=0 . . . N−2, the stored zero-crossing intervals

pulseIntervalMean=Mean(D)

The Shaft Indicators (SI) algorithm 358 calculates miscellaneous shaftindicators included in 358 a. SO1 (Shaft Order 1 in g) is theonce-per-rev energy in the signal average, and is used to detect shaftimbalance. SO2 (Shaft Order 2 in g) is the twice-per-rev energy in thesignal average, and is used to detect shaft misalignment.

SO3 (Shaft Order 3), is the three-per-rev energy in the signal average,and is used to detect shaft misalignment. The miscellaneous shaftindicators may also be included in an embodiment defined as follows:

p=numPathPairs${shaftRatio} = {\frac{\prod\limits_{i = 0}^{p - 1}{shaftPath}_{2i}}{\prod\limits_{i = 0}^{p - 1}{shaftPath}_{{2i} + 1}} = \frac{driving}{driven}}$${indexRatio} = {\frac{\prod\limits_{i = 0}^{p - 1}{indexPath}_{2i}}{\prod\limits_{i = 0}^{p - 1}{indexPath}_{{2i} + 1}} = \frac{driving}{driven}}$${driveRatio} = {\frac{indexRatio}{shaftRatio} \cdot {pulsesUsed}}$${shaftSpeed} = \frac{60}{{pulseIntervalMean} \cdot {driveRatio}}$${resampleRate} = {\frac{shaftSpeed}{60} \cdot {pointsPerRev}}$

RS=residual data,

A1=signal average,${signalAverageL1} = \frac{{P2p}({A1})}{{Rms}\quad ({A1})}$

FF=FFT of the signal average,

shaftOrder_(j)={square root over (FF_(j))}, j=1 . . . 3${gearDistFault} = \frac{{Stdev}\quad ({RS})}{{Stdev}\quad ({A1})}$

As described elsewhere herein, gearDistFault (GDF) is an effectivedetector for distributed gear faults such as wear and multiple toothcracks, and is a complement of the indicator signalAverageL1 (also knownas gearLocalFault).

In addition to the specifically referenced vectors below, the SIalgorithm takes input from the indexer zero-crossing vector (ZC) and mayalso use others and indicated above.

The following definitions for indicators may also be included in anembodiment in connection with the SI algorithm:

shaftPath is defined for the shaft descriptor

indexPath is the path of the shaft seen by the indexer used for signalaveraging

numPathPairs is the number of path pairs defined for shaftPath andindexPath

pulsesUsed is the number of pulses used per revolution of the indexershaft

pulseIntervalMean is the mean of the zero-crossing (ZC) intervals

pointsPerRev is the number of output points per revolution in the signalaverage,

The Bearing Energy (BE) algorithm 376 performs an analysis to reveal thefour bearing defect frequencies (cage, ball spin, outer race, and innerrace frequencies) that usually modulate the bearing shaft frequency. Assuch, these four frequencies are calculated based on the measured shaftspeed and bearing geometry. Alternatively, the four frequency ratios maybe obtained from the bearing manufacturers. The energy levels associatedwith these four frequencies and their harmonics are calculated forbearing fault detection. They are:

Cage Energy: the total energy associated with the bearing cage defectfrequency and its harmonics. Usually it is detectable only at the laterstage of a bearing failure, but some studies show that this indicatormay increase before the others.

Ball Energy: the total energy associated with the bearing ball spindefect frequency and its harmonics.

Outer Race Energy: the total energy associated with the bearing outerrace defect frequency and its harmonics.

Inner Race Energy: the total energy associated with the bearing innerrace defect frequency and its harmonics.

The Total Energy indicator gives an overall measure of the bearingdefect energies.

In one embodiment, one or more algorithms may be used in determining aCI representing a score quantifying a difference between observed oractual test distribution data and a normal probability distributionfunction (PDF) or a normal cumulative distribution function (CDF). Theseone or more algorithms may be categorized as belonging to a class ofalgorithms producing CIs using hypothesis tests (“hypothesis testingalgorithms”) that provide a measure of difference in determining whethera given distribution is not normally distributed. These hypothesistesting algorithms produce a score that is used as a CI. The score maybe described as a sum of differences between an observed or actual testdistribution function based on observed data and a normal PDF or normalCDF. An algorithm may exist, for example, based on each of the followingtests: Chi-Squared Goodness of fit (CS), Kolmogorov-Smimov Goodness offit (KS), Lilliefors test of normality, and Jarque-Bera test ofnormality (JB). Other embodiments may also include other algorithmsbased on other tests for normality, as known to those of ordinary skillin the art. The hypothesis tests compare the test distribution to thenormal PDF, for example as with CS test, or the normal CDF, for exampleas with the KS and Lilliefor tests.

What will now be described is an example in which the CS test is used indetermining a score with a test distribution of observed actual data. Inthis example, the test distribution of observed data forms a Gamma (5,20) distribution function, having and alpha value of 5 and a beta valueof 20. The mean of this Gamma (5,20) distribution is alpha*beta having avariance of alpha*beta². The Gamma (5,20) distribution function is atailed distribution which graphically is similar to that of a normaldistribution.

Referring now to FIG. 19, shown is an example of a graphicalrepresentation 400 of observed data.

Referring now to FIG. 20, shown is an example of a graphicalrepresentation 410 of the normal CDF and the Gamma (5,20) CDF of randomdata. Referring now FIG. 21, shown is an example of a graphicalrepresentation 420 of the difference between the normal CDF and theGamma (5,20) CDF.

In one embodiment, if there are 1000 test samples used in forming asingle CDF, the graphical representation, for example, in FIG. 21represents differences in 1000 instances where the difference betweenthe expected value (Normal CDF) and the maximum deviation of the (inthis case defined as the score) observed gamma CDF can exceed somecritical value. The critical value is that statistic which representssome predefined alpha error (the probability that the test indicates thedistribution is not normal when in fact it is normal—this is typicallyset at 5%.) If the score exceeds the critical value, the distribution issaid to be not normal statistic. The score is the maximum deviation fromthis statistic or alpha value.

It should be noted that the sensitivity or goodness of the testincreases as the number of samples or instances (degrees of freedom “n”)increases approximately as the square root of “n”. For example, in thecase where 1000 instances or samples are used such that n=1000, thesensitivity or ability of this CI to be used in detecting gear faults,for example, is roughly 31 times more powerful than kurtosis inidentifying a non normal distribution.

As another example, in the algorithm using the CS test, the normal PDFis used. Referring now to FIG. 22, shown is a graphical representation430 of the normal PDF and the PDF of the Gamma (5,20) distribution. Therepresentations of FIG. 22 are drawn as continuous lines rather thandiscrete intervals.

Referring now to FIG. 23, the quantities of the x-axis represented inFIG. 22 are shown in another representation 440 as being divided intodiscrete bins, intervals, or categories. For example, there may be 4bins or intervals between any two integer quantities. Between 0 and 1,bin 1 includes values between [0, 0.25), bin 2 includes values between[0.25,0.50), bin 3 includes values between [0.050, 0.75) and bin 4includes values between [0.75,1.0). For each bin, determine the numberof observed and expected values, and their difference. Square each ofthe differences for each bin and then add all the differences and divideby the expected value for each bin. The CS test which sums all thedifferences for each category divided by the expected value for eachcategory represented as:$\sum\limits_{i = 1}^{k}\frac{\left( {{fi} - {ei}} \right)^{2}}{ei}$

for k categories or bins, k−1 degrees of freedom, fi is observed dataand ei is expected data value or number in accordance with a normaldistribution.

For each bin, take the difference between the observed and expectedobservation. Square this value and divided by expected number ofobservation. Sum over all bins. The statistic, the critical value is theχ² at k−1 degrees of freedom may be, for example, 90.72 which is muchgreater than the 0.05 alpha value of a χ², which is 54.57 for 39 degreesof freedom or 40 categories/bins. Thus, the observed data in thisexample as indicated by the statistic is not normally distributed. FIG.24A represents graphically a difference between observed and expectedvalues for each bin or interval of FIG. 23.

It should be noted that the foregoing algorithms provide a way ofmeasuring both the skewness and kurtosis simultaneously by comparing thePDF or CDF of the test distribution against the PDF/CDF of a standardnormal distribution in which a score is used as a CI as described above.

As known to those of ordinary skill in the art, other algorithmsbelonging to the hypothesis testing class may be used in computing CIs.The particular examples, algorithms, and tests selected for discussionherein are representations of those that may be included in the generalclass.

What will now be described is another algorithm that may be used indetermining a CI in an embodiment of the system of FIG. 1. This may bereferred to as an impulse determination algorithm that produces a CIindicating an amount of vibration that may be used in detecting a typeof fault. The impulse determination algorithm takes into account thephysical model of the system. One type of fault that this technique maybe used to detect is a pit or spall on either: gear tooth, inner bearingrace, outer bearing race or bearing roller element. This technique usesa model designed to detect this type of fault where the model is basedon knowledge of the physical system. For example, if there is a pit orspall on a bearing, this may produce a vibration on a first bearingwhich may further add vibrations to other components connected to orcoupled to the bearing.

In one embodiment, a model can be determined for a particularconfiguration by using configuration data, for example. In oneconfiguration, for example, a signal received at a sensor may be asuperposition of gear and bearing noise that may be represented as aconvolution of gear/bearing noise and a convolution of the Gear/Bearingsignal with the gearbox transfer function. Given this, if one type offault is a pit or spall on either a: gear tooth, inner bearing race,outer bearing race or bearing roller element, a model that is designedto look for this type of fault can take advantage of knowledge of thephysical system.

The impulse determination algorithm uses Linear Predictive Coding (LPC)techniques. As known to those skilled in the art, LPC may becharacterized as an adaptive type of signal processing algorithm used todeconvolute a signal into its base components. In the case of apit/spall fault, the base signal components are an impulse traingenerated by the fault hitting a surface (e.g gear tooth with geartooth,inner race with roller element, etc) and the bearing/case transferfunction. The bearing, gear and case have there own transfer functions.Convolution here is transitive and multiplicative. As such, LPCtechniques may be used to estimate the total convolution function of thetotal vibration that may be produced.

For example, in this arrangement, the total amount of vibrationrepresenting the total impulse signal generated by a configuration maybe represented as:

[impulse]f(Gear)f(Bearing)f(Case)≡[impulse][f(Gear)f(Bearing)f(Case)]

in which represents the convolution operation.

It should also be noted that convolution is a homomorphic system suchthat it is monotonically increasing and that logarithmic transformationshold. Thus the relationship of c=a*b also holds for Log c=Log a+Log b. A“dual nature” of convolution is used in following representations toequate operations using convolution in the time domain to equivalentmultiplication operation in the frequency domain.

If “y” represents the total response of all elementary responses, and“h” represents the response of the system for a series of elementaryinput impulses “imp” such that y is the convolution of imp and h, thenthis may be represented as:

y=imph

and then converting “y” and “h” each, respectively, to the frequencydomain represented as “Y” and “H”, as may be represented by thefollowing:

Y=ℑ(y), H=ℑ(h)

taking the Fourier transform (FFT) of each where H represents thetransfer function. The convolution in the time domain may be equated toa multiplication in the frequency domain represented as:

Y=IMP·H

in which IMP is the Fourier transformation of imp into the frequencydomain. Above, imp is in the time domain.

The convolution in the time domain is equivalent to multiplication inthe Frequency Domain. Referring to the homomorphic property ofconvolution, it follows that:

 log(Y)=log(IMP)+log(H),

therefore

log(IMP)=log(Y)−log(H),

IMP=exp(log(Y)−log(H))

and finally

imp=ℑ⁻¹(IMP)

Using the foregoing, the system transfer function “H” may be estimatedfor the Gear/Bearing and Case to recover the impulse response allocatedwith a Gear or Bearing pit/spall fault. The estimation of this transferfunction “H” may be accomplished using Linear Predictive Coding (LPC)techniques. LPC assumes that the Transfer Function is a FIR filter, andas such, the auto-correlation of the time domain signal may be used tosolve for the filter coefficients in a minimum sum of square errorsense.

Using the LPC model, there is an impulse that is convoluted with a FIRfilter, such that:

y[n]=a ₁ x[n−1]+a ₂ x[n−2]+a ₃ x[n−3]+ . . .

LPC techniques may be used to estimate the coefficients a=(a₁ . . . an)for an order p in a minimum sum of square error sense, n=p+1. Thestandard least squares error estimators may be used, wherein y=y[1, 2, .. . n], and x is the time delayed signal, in which: $x = \begin{bmatrix}\begin{matrix}{x\left\lbrack {{n - 1},{n - 2},{{\ldots \quad n} - p}} \right\rbrack} \\{x\left\lbrack {{n - 2},{n - 3},{n - p - 1}} \right\rbrack}\end{matrix} \\\begin{matrix}\vdots & \vdots\end{matrix}\end{bmatrix}$

where a=(x^(T)x)⁻¹x^(T)y. These values for a1 . . . an may be used withthe following equation:

y_(hat)=ax, b=(y−y_(hat))² and the estimator of error B is: Σ_(all)b.

Y may also be expressed as:

Y=FFT(y[1, 2, . . . n])

in which y[1 . . . n] are values in the time domain expressed in thefrequency domain as a Fourier transform of the time domain values. Yrepresents current time vector measurements in the frequency domain.

In terms of a and B, the transfer function H may be estimated andrepresented as a/B, (freq. Domain). Note that “a” is a vector of thevalues a1 . . . an obtained above.

The homomorphic property of convolution as described above may be usedto estimate the impulse as represented in:

IMP=exp(log(Y)−log(H)) IMP Equation

If there is no fault, the impulse, for example, may be characterized as“white noise”. As the fault progresses, the impulse or the value of Hbecomes larger. The CI is the power spectral density at a bearingpassing frequency for a bearing fault, or a mesh frequency for a gearfault. Other CIs based on the foregoing value may be a “score” of theLilifers test for normality, or other such test.

In the foregoing, a pit or spall may cause a vibration or tapping.Subsequently, other elements in contact with the ball bearing may alsovibrate exhibiting behavior from this initial vibration. Thus, theinitial vibration of the pit or spall may cause an impulse spectrum tobe exhibited by such a component having unusual noise or vibration.

The value of IMP as may be determined using the IMP Equation aboverepresents the impulse function that may be used as a “raw” value and ata given frequency and used as an input into an HI determinationtechnique. For example, the IMP at a particular frequency, since thisthe spectrum, determined above may be compared to expected values, suchas may be obtained from the stored historic data and configuration data.An embodiment may also take the power spectrum of this raw impulsespectrum prior to being used, for example, as input to an HI calculationwhere the power spectrum is observed at frequencies of interest, such asthe inner race frequency. For example, if the impulse function is withinsome predetermined threshold amount, it may be concluded that there isno fault.

What is shown in the FIG. 24B and FIG. 24C are relative to a healthysystem, such as a main gearbox, for example, such as in connection witha planetary race fault of an SH-60B U.S. Navy Helicopter built bySilorsky.

The FIG. 24B representation 700 shows an impulse train in the frequencydomain of the healthy system.

It should be noted that an embodiment may estimate the transfer functionH using LPC using different techniques. An embodiment may estimate thetransfer function H using an autocorrelation technique (AutoLPC). Anembodiment may also estimate the transfer function H using a covariancetechnique (CovLPC). Use of autocorrelation may use less mathematicaloperations, but require more data than using the covariance.Alternatively, use of the covariance technique may use more mathematicaloperations but require less data. As the amount of available dataincreases, the autocorrelation LPC result converges to the covarianceLPC result. In one example, data samples are at 100 KHz with 64,000 datapoints used with the autocorrelation technique due to the relativelylarge number of data points.

FIG. 24C representation 710 shows the data of 700 from FIG. 24B in thetime domain rather than the frequency domain.

FIG. 24D representation 720 shows the power spectral density of theabove figures as deconvolved time data of frequency v. dB values in ahealthy system.

The foregoing FIGS. 24B-24D represent data in a graphical display inconnection with a healthy system. Following are three additionalgraphical displays shown in FIGS. 24E-24G in connection with anunhealthy system, such as a starboard ring channel which exhibit datathat may be expected in connection with a pit or spall fault.

FIG. 24E, representation 730, illustrates an impulse train as may beassociated with an unhealthy system in the time domain. FIG. 24F,representation 740, illustrates a graphical display of the impulse trainin the frequency domain.

In FIG. 24G, shown is an illustration 740 is a graphical representationof the power spectrum of the impulse train represented in connectionwith the other two figures for the unhealthy system identified by aperiod impulse train associated with an inner race bearing fault. Inthis example, a spike may be viewed in the graphical display as well asthe harmonics thereof.

It should be noted that other algorithms and CIs in addition to thosedescribed herein may be used in producing CIs used in techniques inconnection with HIs elsewhere herein.

What will now be described is one embodiment in which these CIs may beused. Referring now to FIG. 25, shown is a flow chart of steps of oneembodiment for determining the health of a part as indicated by an HI.At step 502, raw data acquisition is performed. This may be, forexample, issuing appropriate commands causing the VPU to perform a dataacquisition. At step 504, the raw data may be adjusted, for example, inaccordance with particular configuration information producing analysisdata as output. It is at step 504, for example, that an embodiment maymake adjustments to a raw data item acquired as may be related to theparticular arrangement of components. At step 506, data transformationsmay be performed using the analysis data and other data, such as rawdata. The output of the data transformations includes transformationoutput vectors. At step 508, CIs are computed using the analysis dataand transformation vector data as may be specified in accordance witheach algorithm. At step 510, one or more CIs may be selected. Particulartechniques that may be included in an embodiment for selectingparticular CIs is described elsewhere herein in more detail. At step512, CIs may be normalized. This step is described in more detailelsewhere herein. At step 514, the selected and normalized CIs are usedin determining HIs. Particular techniques for determining HIs aredescribed in more detail elsewhere herein.

In an embodiment, due to the lengthy processing times, for example, inexecuting the different algorithms described herein, HI computations maynot be executed in real time. Rather, they may be performed, forexample, when a command or request is issued, such as from a pilot or atpredetermined time intervals.

The hardware and/or software included in each embodiment may vary in oneembodiment, data acquisition and/or computations may be performed by oneor more digital signal processors (DSPs) running at a particular clockspeed, such as 40 MHz, having a predetermined numerical precision, suchas 32 bits. The processors may have access to shared memory. In oneembodiment, sensors may be multiplexed and data may be acquired ingroups, such as 8. Other embodiments may vary the number in each groupfor data sampling. The sampling rates and durations within anacquisition group may also vary in an embodiment. Data may be placed inthe memory accessed by the DSPs on acquisition. In one embodiment, thesoftware may be a combination of ADA95 and machine code. Processors mayinclude the VPU as described herein as well as a DSP chip.

What will now be described are techniques for normalizing CIs inconnection with determining HIs providing more detailed processing ofstep 512 as described in connection with flowchart 500.

Transmission error (T.E.) depends upon torque. Additionally, vibrationdepends upon the frequency response of a gear. As such, the CI, whichalso depends upon T.E. and vibration, is a function (generally linear)of torque and rotor speed (which is frequency), and airspeed as this maychange the shape of the airframe. Thus, techniques that may be used inconnection with determining the “health state” or HI of a component maynormalize CIs to account for the foregoing since HIs are determinedusing CIs.

For each bearing, shaft and gear within a power train, a number of CIsmay be determined. An embodiment may compare CI values to thresholdvalues, apply a weighting factor, and sum the weighted CIs to determinean HI value for a component at a particular time. Because dataacquisitions may be made at different torque (e.g. power setting)values, the threshold values may be different for each torque value. Forexample, an embodiment may use 4 torque bands, requiring 4 thresholdvalues and weights for each CI. Additionally, the coarseness of thetorque bands will result in increased, uncontrolled system variance.Alternatively, rather than use multiple threshold values and have anuncontrolled variance, an embodiment may use a normalization techniquewhich normalizes the CI for torque and rotor RPM (Nr), and airspeed,expressed as a percentage, for example, in which a percentage of 100% isperfect. Use of these normalized CIs allows for a reduction ofconfiguration such that, for example, only one threshold is used andvariance may also be reduced.

The normalization technique that will now be described in more detailmay be used in connection with methods of HI generation, such as thenon-linear mapping method and the hypothesis testing method of HIgeneration that are also described in more detail elsewhere herein.

It should be noted that a deflection in a spring is linearly related tothe force applied to the spring. The transmission may be similar incertain aspects to a large, complex spring. The displacement of a pinionand its corresponding Transmission Error (T.E.) is proportional to thetorque applied. T.E. is a what causes vibration, while the intensity ofthe vibration is a function of the frequency response (N_(r)), wherefrequency is a function of RPM. Thus, vibration and the corresponding CIcalculated using a data acquisition are approximately linearlyproportional to torque, N_(r). (over the operating range of interest)and/or airspeed although at times there may be a linear torque*Nrinteraction effect. For example, gear box manufacturers may design agearbox to have minimum T.E. under load, and a graphical representationof T.E. vs. Torque is linear, or at least piece wise linear. It shouldbe noted that test data, for example used in connection with a Bellhelicopter H−1 loss of lube test, shows a relationship between CI andtorque suggesting linearity. Additionally, tests show that airspeed isalso relevant factor. Other embodiments may take into account any one ormore of these factors as well as apply the techniques described hereinto other factors that may be relevant in a particular embodiment orother application although in this example, the factors of torque,airspeed and Nr are taken into account.

An equation representing a model minimizing the sum of square error of ameasured CI for a given torque value in a healthy gear box is:

CI=B ₀ +B ₁*Torque+B ₂ Nr+B ₃Airspeed+T.E  (Equation 1)

The order of the model may be determined by statistical significance ofthe coefficients of Equation 1. In the previous equation, the T.E. of a“healthy” component may have, for example, a mean of zero (0) with someexpected variance. It should be noted that if the model fits well forthe lower order. Higher order coefficients are not required and mayactually induce error in some instances. The following example is builtas a first order model, higher orders may be solved by extension of thatexplained in the first order model. This model, written in matrix formatis: y=B x where $y = {{\begin{bmatrix}{CI}_{1} \\{CI}_{\ldots} \\{CI}_{n}\end{bmatrix}\quad B} = {{\begin{bmatrix}B_{0} & {B_{1}\quad \ldots \quad B_{n}}\end{bmatrix}\quad \text{and}\quad x} = \begin{bmatrix}1 & t_{1} & N_{R_{1}} & {Airspeed1} \\1 & {t\quad \ldots} & N_{R\quad \ldots} & {{Airspeed}\quad \ldots} \\1 & t_{n} & {N_{R}n} & {Airspeedn}\end{bmatrix}}}$

Each of the CIs included in the vector y is a particular recorded valuefor a CI from previous data acquisitions, for example, as may be storedand retrieved from the collected data 18. Also stored with eachoccurrence of a CI for a data acquisition in an embodiment may be acorresponding value for torque (t), Nr, and Airspeed. These values mayalso be stored in the collected data 18.

The model coefficients for B may be estimated by minimizing the sum ofsquare error between the measured CI and the model or estimated CI usingthe observed performance data. Solving the foregoing for the unbiasedestimator of B=(x^(T)x)⁻¹x^(T)y. The variance of B is:Var(B)=E(b−B)(b−B)^(T)=σ²(x^(T)x)⁻¹ where b is an unbiased estimator ofB. The unbiased estimator of σ² is s²:$s^{2} = {\frac{e^{T}e}{n - p - 1} = {\frac{\left( {y - \hat{y}} \right)^{T}\left( {y - \hat{y}} \right)}{n - p - 1} = \frac{{y^{T}y} - {b^{T}x^{T}y}}{n - p - 1}}}$

In the vector B from y=xB, coefficient B₀ represents the mean of thedata set for a particular component which, for example, may berepresented as an offset value. Each of the other values B1 . . . Bn arecoefficients multiplied by the corresponding factors, such as airspeed,torque, and Nr.

The foregoing B values or coefficients may be determined at a time otherthan in real-time, for example, when flying a plane, and thensubsequently stored, along with corresponding X information, forexample, in the collected data store 18. These stored values may be usedin determining a normalized CI value for a particular observed instanceof a CIobs in determining an HI. The normalized CI may be representedas:

CI _(normalized) =T.E.=CIobs−(B*x)

where CIobs represents an instance of a CI being normalized usingpreviously determined and stored B and x values. Threshold values, asmay be used, for example, in HI determination, may be expressed in termsof multiples of the standard deviation Warning=B₀+3*σ²(x^(t)x)⁻¹,Alarm=B₀+6*σ²(x^(t)x)⁻¹. It should be noted that a covariance that maybe determined as:

Σ=s²(x^(t)x)⁻¹ where s² is calculated as noted above.

As described elsewhere herein, the foregoing techniques are based upon ahealthy gear characterized as having noise that is stationary andGaussian in which the noise approximates a normal distribution.

What will now be described are techniques that may be used indetermining an HI using the normalized CI values as inputs. Inparticular, two techniques will be described for determining an HI. Afirst technique may be referred to as the non-linear map technique. Thesecond technique may be referred to as the hypothesis test method of HIgeneration. It should be noted that CI values other than normalized CIvalues may be used in connection with HI determination techniquesdescribed herein.

It should be noted that an embodiment may use CI values that are notnormalized in connection with the HI determination techniques describedherein. In this instance, multiple torque bands may be used, one foreach CI or group of CIs belonging to different torque bands.Additionally, a larger covariance matrix may be used as there may be alarger variance causing decrease in separation between classes.

For any generic type of analysis (gear, bearing, or shaft), a subset ofthe diagnostics indicators or CIs is selected. The CIs which are bestsuited to specify the fault indication may be developed over timethrough data analysis. Faults may be calculated at the component leveland an HI may be calculated for a given component. If there is acomponent fault, then there is a sub-assembly fault, and therefore adrive train fault.

Following is a description of a non-linear mapping methodology fordetermining an HI. Given a set of component indicators I1, I2, I3, . . .IN, choose the desired subset of K indicators such that K<=N. For thechosen group of indicators, let WTi define the weight of the ithindicator, Wi the warning threshold, and Ai the alarm threshold. Thenapply the following processing to the set of chosen indicators.

Health Indicator Contribution Description

for XX = 1:K /* cycle through all K indicators in subset */ If I[XX] <Wi /* if less than warning level Wi, assign 0 */ Hi contribution = 0elseif Wi * Ii < Ai Hi contribution = 1*Wi else Hi contribution = 2*Wiend end

In the foregoing pseudo-code like description, each indicator or CI isweighted and contributes a portion to the HI determination. Subsequentlyall the Hi contributions for the selected CIs are summed and may becompared to threshold values for determining one of two possibleoutcomes of “healthy” or “not healthy”.

Consider the following example table of information for a selectedsubset of 9 CIs along with threshold and weight values. It should benoted that in an embodiment, any one or more of the values for weights,warning and alarm values may be modified.

CI Warning Alarm No. Value Level Level Weight HI contribution I2 3.263.5 4.0 1.0 0.0 I3 3.45 3.0 3.5 1.0 1.0 I6 7.5 6.0 8.0 1.4 1.4 I9 0.880.5 0.75 0.9 1.8 I14 4.2 3.5 4.5 1.0 1.0 I17 4.7 3.5 4.5 0.9 1.8 I22 5.22.0 4.0 1.1 2.2 I23 4.4 3.5 4.5 1.2 1.2 I24 18.9 10.0 20.0 1.0 1.0

Using the foregoing example and values, the sum of the HI contributionsis 11.4. Applying the Health Indicator Contribution technique as setforth in the foregoing pseudo-code like description, I2, with a value of3.26, is below the warning threshold, so the contribution to the indexis 0. Indicator I3 has a value of 3.45, which contributes a 1 toward theindex since the weight value is also 1. However, Indicator I6contributes a 1.4 to the index because it crosses the warning level(contributing a value of 1 to the index) while being weighted by afactor of 1.4.

In the foregoing example, if no indicators were in alarm, the sum of HIcontributions would be zero and if all indicators were in alarm, the sumwould be 19, the worst fault case represented by this detector scheme.The HI may be represented as a value of 1 for healthy and 0 for nothealthy as associated with a component represented by the foregoing CIvalues.

The HI may be determined by dividing 11.4/19, the maximum of worst caseoutcome to obtain 0.6. This overall health index output ratio can thenbe compared to another final output threshold, where normal componentsproduce His, for example, less than 0.5; values between 0.5 and 0.75represent warning levels, and values over 0.75 represent alarm.

It should be noted that the weights may be determined using a variety ofdifferent techniques. The weights of each CI may be determined using anyone or more of a variety of techniques. One embodiment may determineweights for the CIs as:$\frac{1}{\sqrt{{eigen\_ values}{\_ of}{\_ the}{\_ covariance}{\_ matrix}}}$

It should be noted that other threshold values may be used in HIdetermination and may vary with each embodiment.

In one embodiment, using the normalized CI described elsewhere hereinwith the non-linear mapping technique, the threshold values may berepresented as: Warning=B₀+3*σ²(x^(t)x)⁻¹, Alarm=B₀+6*σ²(x^(t)x)⁻¹,where B₀ may represent a mean or average coefficient as included in theB vector being solved for in the equations described in connection withCI normalization. In the foregoing example, the Warning threshold is 3standard deviations and the Alarm level is 6 standard deviations. Itshould be noted that other threshold values may be used in and may varyin accordance with each embodiment.

What will now be described is a second technique that may be used indetermining HIs using CIs, in particular, using normalized CIs.

The technique for HI determination may be referred to as Hypothesistesting technique for HI determination which minimizes the occurrence ofa false alarm rate, or incorrectly diagnosing the health of a part asbeing included in the alarm classification when in fact the part is notin this particular state. In one embodiment, three classes of healthindication may be used, for example, normal, warning and alarmclassifications with alarm being the least “healthy” classification.Other embodiments may use the techniques described herein with adifferent number of classes. As described elsewhere herein, the class ofa part indicating the health of the part may be determined based onmeasured vibrations associated with the part. Additionally, thetechnique described herein may use a transformation, such as thewhitening transformation to maximize the class distributions orseparation of values thus decreasing the likelihood or amount of overlapbetween the classes. In particular, this maximization of classseparation or distance attempts to minimize the misclassification of apart. A description of the whitening transformation used in herein infollowing paragraphs may be found, for example, in “Detection,Estimation and Modulation Theory”, Harry L. Van Trees, 1968, John Wiley& Sons, New York Library of Congress Catalog Card Number 67-23331.

Using the Hypothesis Testing method of HI generation, the HI orclassification h(X) of a vector of normalized CI values denoted as X maybe determined in which, as discussed elsewhere herein in more detail, Xmay be normalized. Using the hypothesis testing technique, adetermination is made as to which class (normal, warning or alarm) Xbelongs. In our instance, there are three classes. However, a firstdetermination using the hypothesis testing may be performed using afirst class corresponding to normal, and a second class corresponding tonot normal. If the determination is normal, then testing may stop.Otherwise, if determination is made that the testing results are “notnormal”, a further or second determination using the hypothesis testingmay be performed to determine which “not normal” class (alarm orwarning) X belongs. Thus, the hypothesis testing technique may beperformed more than once in accordance with the particular number ofclasses of an embodiment. For three classes, there are two degrees offreedom such that if the sample X is not from A or B classes, then it isfrom Class C.

X may belong to class ω₁ or ω₂, such that:${q_{1}(X)}\quad \underset{\omega_{2}}{\overset{\omega_{1}}{\Diamond}}\quad {q_{2}(X)}$

(the notation $\underset{\omega_{2}}{\overset{\omega_{1}}{\Diamond}}$

means that if q₁(X) is greater than q₂(X), choose class 2, ω₂, or ifq₁(X) is less than q₂(X), choose class 1, ω₁.) In the foregoing, q_(i)is the a posteriori probability of ω_(i) given X, which can be computed,using Bayes theorem in which q_(i)=P_(i)p_(i)(X)/p(X), where p(X) is themixed density function.

The mixed density function is the probability function for all caseswhere q_(i) is the unconditional probability of “i” given theprobability of “i” conditioned on the mixed density function.

Substituting the foregoing representation of each q1 and q2, since p(X)is common to both, now:$P_{1}{p_{1}(X)}\quad \overset{\omega_{1}}{\underset{\omega_{2}}{\Diamond}}\quad P_{2}{p_{2}(X)}$

or as a likelihood function as${l(X)} = {\frac{p_{1}(X)}{p_{2}(X)}\quad \overset{\omega_{1}}{\underset{\omega_{2}}{\Diamond}}\quad {\frac{P_{2}}{P_{1}}.}}$

The likelihood ratio is a quantity in hypothesis test. The value P₂/P₁is the threshold value. In some instances, it may be easier to calculatethe minus log likelihood ratio. In this case, the decision rule becomes(e.g. now called the discriminate function):${h(X)} = {{{- \ln}\quad {l(X)}} = {{{- \ln}\quad {p_{1}(X)}} + {\ln \quad {p_{2}(X)}\quad \overset{\omega_{1}}{\underset{\omega_{2}}{\Diamond}}\quad \ln \quad \frac{P_{2}}{P_{1}}}}}$

Assume that the p_(i)(X)'s are normally distributed with mean orexpected values in vectors M_(i′) and covariance matrix Σ_(i). Thisassumption may be determined without loss of generality in that, anynon-normal distribution can be whitened, as with the whiteningtransformation described elsewhere herein, with the appropriate powertransform, or by increasing the sample size to the point where thesample size is very large. Given this, the decision rule becomes:$\begin{matrix}\begin{matrix}{{h(X)} = \quad {{- \ln}\quad {l(X)}}} \\{= \quad {{\frac{1}{2}\left( {X - M_{1}} \right)^{T}{\overset{- 1}{\sum\limits_{1}}\left( {X - M_{1}} \right)}} -}} \\{\quad {{\frac{1}{2}\left( {X - M_{2}} \right)^{T}{\sum\limits_{2}^{- 1}\left( {X - M_{2}} \right)}} +}} \\{\quad {\frac{1}{2}\ln \quad \frac{\sum\limits_{1}}{\sum\limits_{2}}\quad \overset{\omega_{1}}{\underset{\omega_{2}}{\Diamond}}\quad \ln \frac{P_{2}}{P_{1}}}}\end{matrix} & \text{Equation~~E1}\end{matrix}$

Recall that maximization of distance between the two classes is desiredto minimize the chance of a false alarm or misclassification of a partas broken when it is actually normal.

A function Z is defined as Z=X−M, (e.g. a shift where X is the measuredCI data and M is the mean CI values for a class), so that:

d_(z) ²(z)=Z^(T)Σ⁻¹Z (this distance is the n dimensional distancebetween two distributions). Note that Σ represents the covariance. Itmay be determined that a particular Z maximizes the distance function,subject to Z^(T)Z=I, the identity matrix.

Using a standard Lagrange multiplier, μ, to find the local extrema (e.g.the maximum) a partial derivative is obtained with respect to Z in thefollowing:${\frac{\partial\quad}{\partial Z}\left\{ {{Z^{T}\Sigma^{- 1}Z} - {\mu \left( {{Z^{T}Z} - I} \right)}} \right\}} = {{2{\sum\limits^{- 1}Z}} - {2\mu \quad Z}}$

where Σ is the eigenvector of X,

which may then be set to zero to find the extrema and solving for Z:Σ⁻¹Z=μZ or ΣZ=λZ where λ={fraction (1/μ)}. In order that a non-null Zexits, λ must be chosen to satisfy the determinant: |Σ−λI|=0.

Note that λ is the eigenvalue of X and Σ is the correspondingeigenvector. Σ is a symmetric n×n matrix (e.g. a covariance matrix),there are n real eigenvalues (λ₁ . . . λ_(n)) and n real eigenvectors φ₁. . . φ_(n). the characteristic equation is: ΣΦ=ΦΛ, and Φ^(T)Φ=I where Φis an n×n matrix consisting of n eigenvectors and Λ is a diagonal matrixof eigenvalues (e.g. the eigenvector matrix and eigenvalue matrix,respectively).

Y, representing the coordinated shifted value of X, may be representedas:

Y=Φ^(T)X,

having a covariance matrix of y, Σ_(y)=Φ^(T)Σ_(x)Φ=Λ where Σ_(x)represents the covariance of the vector of matrix x. Continuing, thewhitening transformation may be defined such that:

 Y=Λ ^(−1/2)Φ^(T) X=(ΦΛ^(−1/2))^(T) X, Σ_(y)=Λ^(−1/2)Φ^(T)Σ_(x)ΦΛ^(−1/2)=Λ^(−1/2)ΛΛ^(−1/2)=I,

Thus the transformation that maximizes that distance betweendistribution or classes is:

A=Λ^(−1/2)Φ^(T)

as shown above.

Using this value of A, define${{A^{T}{\sum\limits_{1}A}} = I},{{A^{T}{\sum\limits_{2}A}} = K},\quad {and}$A^(T)(M₂ − M₁) = L  and

(Σ₁ ⁻¹Σ₂ ⁻¹)⁻¹ transformed to a diagonal matrix Λ by A that may berepresented as:

Λ=A ^(T) [A(I−K ⁻¹)A ^(T)]⁻¹ A=(I−K ⁻¹)⁻¹

which may be substituted into the discriminate function defined above:${h(X)} = {{\frac{1}{2}Y^{T}\Lambda^{- 1}Y} - {\left\lbrack \left( {{- K^{- 1}}L} \right)^{T} \right\rbrack Y} + \left\lbrack {{{- \frac{1}{2}}L^{T}K^{- 1}L} - {\frac{1}{2}\ln \quad {K}} - {\ln \frac{P_{2}}{P_{1}}}} \right\rbrack}$

Thus, if the above is less than the threshold, for example, ln(P₂/P₁),then the component is a member of the normal or healthy class.Otherwise, the component is classified as having an HI in the brokenclass, such as one of alarm or warning. In the latter case, anotheriteration of the hypothesis testing technique described herein may befurther performed to determine which “broken” classification, such asalarm or warning in this instance, characterizes the health of thecomponent under consideration.

In the foregoing technique for hypothesis testing, values, such as the aposteriori probabilities q₁ and q₂, may be obtained and determined priorto executing the hypothesis testing technique on a particular set of CInormalized values represented as X above. As known to those of ordinaryskill in the art, Bayes theorem may be used in determining, for example,how likely a cause is given that an effect has occurred. In thisexample, the effect is the particular CI normalized values and it isbeing determined how likely each particular cause, such as a normal orbroken part, given the particular effects.

It should be noted that operating characteristics of a system define theprobability of a false alarm (PFA) and the probability of detection(PD). The transformation used to maximize the distance functionoptimizes the discrimination between classes. However, the thresholdvalue selected given a discriminate function may be used in determiningthe PD and PFA. In some embodiments, the cost of a false alarm may behigher than the cost of a missed detection. In these instances, the PFAmay be set to define threshold values, and then accept the PD (e.g., aconstant false alarm rate (CFAR) type of process). The distance functionis a normal density function, based on the conditional covariance of thetested values under consideration. Given that, the PFA may be determinedas: P_(F)=P(HoH₁), which means the probability that the sufficientstatistic is greater than some threshold is the integral of thethreshold to infinity of a normal PDF.$P_{FA} = {{\int_{\alpha}^{\infty}{{p_{l|H_{o}}\left( l \middle| H_{o} \right)}{L}}} = {\int_{\alpha}^{\infty}{\frac{1}{2\pi}\exp \quad \left( {- \frac{x^{2}}{2}} \right){x}}}}$

where

the lower integral limit of

α=ln(P ₁ /P ₂)/d+d/2,

and, as before

d ²=(M ₂ −M ₁)^(T)Σ₁ ⁻¹(M ₂ −M ₁)

In this example, the threshold may be the ln(P₂/P₁). This integration isthe incomplete gamma function. Conversely, the probability of adetection (PD) is:$P_{D} = {{\int_{- \infty}^{\alpha}{{p_{l|H_{1}}\left( l \middle| H_{1} \right)}{L}}} = {\int_{- \infty}^{\alpha}{\frac{1}{2\pi}\exp \quad \left( {- \frac{(d)^{2}}{2}} \right){x}}}}$

but now

α=−ln(P ₂ /P ₁)/d+d/2,

and

d ²=(M ₁ −M ₂)^(T)Σ₂ ⁻¹(M ₁ −M ₂)

Note, the distance function is relative to the condition (e.g. H₀ or H₁)being investigated.

Referring now to FIG. 26, shown is an example of a graphicalillustration of the probability of a false alarm PFA represented by theshaded region A3 which designates the overlap between the distributionof class H0, denoted by the curve formed by line A1, and class H1,denoted by the curve formed by line A2.

Referring now to FIG. 27, shown is an example of a graphicalillustration of the probability of an appropriate detection (PD)represented as area A4 as belonging to class represented by H1 asrepresented by the curve formed by line A2.

Referring now to FIG. 28, shown is a graphical illustration of arelationship in one embodiment between the PFA and PD and the thresholdvalue. Note that as the threshold increases, the PD increases, but alsothe PFA increases. If the performance is not acceptable, such as the PFAis too high, an alternative is to increase the dimensionality of theclassifier, such as by increasing the population sample size, n. Sincethe variance is related by 1/sqroot(n), as n increases the variance isdecreased and the normalized distance between the distributions willincrease. This may characterize the performance of the system. Thelikelihood ratio test used herein is a signal to noise ratio such thatthe larger the ratio, (e.g., the larger the distance between the twodistributions), the greater the system performance. The process oftaking an orthonormal transformation may be characterized as similar tothe of a matched filter maximizing the signal to noise ratio.

Referring now to FIG. 29, shown is an example of a graphicalillustration of how the threshold may vary in accordance with theprobability of determining class Ho.

It should be noted that false alarm rate and detection rate are twofactors that may affect selection of particular values, such asthresholds within a particular system. In the example embodimentdescribed herein, false alarm rate is a determining factor, for example,because of the high cost associated with false alarms and the fact thatthey may corrode confidence when a real fault is detected. It should benoted that other embodiments and other applications may have differentconsiderations. Further in this example of the system of FIG. 1, certainfactors may be considered. An acceptable false alarm rate, for example,such as 1 false alarm per 100 flight hours, is established. An estimateof the number of collection opportunities per flight hours may bedetermined, such as four data collections. A number of His may beselected for the system, such as approximately 800. A confidence levelmay be selected, such as that there is a 90% probability that a falsealarm rate is less than 1 per 100 flight hours.

In this example, it should be noted that each HI is a an independentclassification event such that the law of total probability may give thesystem alarm rate using the foregoing:

System PFA=1/(100*4*800)=3.1250*10⁻⁶.

It should also be noted that in the foregoing, when the covariance oftwo classes is approximately the same, or for example, unknown for aclass, the logarithm likelihood ratio test for classification may besimplified in that the model may be reduced to a linear rather thanquadratic problem having the following model:${\left( {M_{2} - M_{1}} \right)^{T}{\sum\limits^{- 1}X}} + {\frac{1}{2}\left( {{M_{1}^{T}{\sum\limits^{- 1}M_{1}}} - {M_{2}^{T}{\sum\limits^{- 1}M_{2}}}} \right)\quad \overset{\omega_{1}}{\underset{\omega_{2}}{\Diamond}}\quad \ln \quad \frac{P_{2}}{P_{1}}}$

If the covariance is whitened, the model simplifies further (assumingthe appropriate transformation is made to the means and measuredvalues).${\left( {M_{2} - M_{1}} \right)^{T}X} + {\frac{1}{2}\left( {{M_{1}^{T}M_{1}} - {M_{2}^{T}M_{2}}} \right)\quad \overset{\omega_{1}}{\underset{\omega_{2}}{\Diamond}}\quad \ln \quad \frac{P_{2}}{P_{1}}}$

What will now be described are techniques that may be used in connectionwith selecting a subset of CIs, such as selection of normalized CIs, forexample, under consideration for use in determining a particular HI.

If we have a two or more classes (such as alarm, warning and normalclassifications), feature extraction, or determining which CIs to use inthis embodiment, may become a problem of picking those CIs or featuresthat maximize class separability. Note that separability is not adistance. As described elsewhere herein, an eigenvector matrixtransformation may be used in maximizing the distance between twofunctions or distribution classes. However, this same technique may notbe applicable when some of the information (e.g. dimensionality) isbeing reduced. For example, in the following test case, three features,or CIs, are available, but only two are to be selected and used indetermining HI classification. The distributions are:${{Cov}_{1} = \begin{bmatrix}1 & {- {.5}} & {.5} \\{- {.5}} & 2 & {.8} \\{.5} & {.8} & 2.5\end{bmatrix}},{M_{1} = \begin{bmatrix}0 \\0 \\0\end{bmatrix}},{{Cov}_{2} = \begin{bmatrix}1 & {.7} & {.7} \\{.7} & 2.5 & 1 \\{.7} & 1 & 2.5\end{bmatrix}},{M_{2} = \begin{bmatrix}3 \\{- 1} \\3\end{bmatrix}}$

When looking at the eigenvalues of the whitening transformation (1.9311,3.0945, 0.4744), the maximum distance of the distribution is an axis y(e.g. 2^(nd) dimension, the distribution was whitened and the projectdimension (e.g. x, y or z) was plotted), but this axis has the minimumseparability. Using this as one of the two features will result inhigher false alarm rates than another feature. This may identify theimportance of feature selection in maximizing the separability.

The problem of separability may be characterized as a “mixed” problem inthat differences in means may be normalized by different classcovariance. If the mean values are the same, or the covariance are thesame, techniques such as the Bhattacharyya Distance may be used tomeasure class separability. However, same mean or covariance values maynot be likely and thus such techniques may not be applicable.Statistical tools developed in discriminant analysis may be used toestimate class separability.

A measure of within class scatter may be represented as the weightedaverage of the class covariance:${S_{w} = {\sum\limits_{i = 1}^{L}{P_{i}\sum\limits_{i}}}},$

for each class I, where Pi is the probability of the occurrence of thecovariance Σ_(I) for that class. In one embodiment, there may be twoclasses, such as healthy or unhealthy. When considering the unhealthystatus, for example, when performing a second round of hypothesistesting described herein, there may be alarm and warning classes.

A measure of between class scatter, Sb, may be represented as themixture of class means:${S_{b} = {\sum\limits_{i = 1}^{L}{{P_{i}\left( {M_{i} - M_{0}} \right)}\left( {M_{i} - M_{0}} \right)^{T}}}},{M_{0} = {\sum\limits_{i = 1}^{L}{P_{i}{M_{i}.}}}}$

Note that M₀ represents the mean or expected value of the classes andMi−M₀ is a difference or variation from the expected value for theclasses under consideration. The formulation for a criteria for classseparability may result in values that are larger when the between classscatter is larger, or when the within class scatter is smaller. Atypical criteria for this is J=diag(S_(w) ⁻¹S_(b)), where In general,S_(w) is not diagonal. One technique takes the whitening transformationof S_(w) where A^(T) S_(w)A=I, then define the whitening transformationof Sb as:

S _(bw) =A ^(T) S _(b) A.

Now taking the diagonal of the foregoing Sbw gives a betterrepresentation of the class separability of each feature.

In summary, CIs may be selected in accordance with the techniquedescribed above to obtain and examine the diagonals of the “whitened”Sb, represented as Sbw. Let X be a matrix where rows and columnsrepresent different CIs having a covariance matrix Σ. An embodiment mayuse normalized CIs and select a portion of these for use. An embodimentmay also use CIs however, those selected should belong to the sametorque band.

As described elsewhere herein, let A represent the correspondingeigenvalue matrix and Φ as the corresponding eigenvector matrix for theCI matrix X. Then, A, as described elsewhere herein in connection withthe whitening transformation, may be represented as:

 A=Λ ^(−1/2)Φ^(T)

where A is the transformation matrix that whitens the covariance Σ. IfSb is defined as above as the between mean covariance of the classes,the whitening matrix A may be used to normalize the differences and givea distance between the mean values of the different classes, such that

Sbw=A ^(T) SbA

where Sbw represents the “whitened” Sb. The diagonals of Sbw may then besorted in descending order in which each diagonal represents anapproximation of the size of the separation between features or CIs.Thus, selection of a subset of “n” features or CIs from a possible “m”maximum CIs included in X may be determined by selecting the “n” largestdiagonals of the matrix Sbw. In particular, the diagonal entry 1,1corresponds to the first column of the covariance matrix and the firstCI in the vector X, entry 2,2 to the second column of the covariancematrix and the second CI in the vector X being considered, and so on.

Once a particular HI is determined at a point in time, it may be desiredto use techniques in connection with trending or predicting HI values ofthe component at future points in time. Techniques, such as trending,may be used in establishing, for example, when maintenance orreplacement of a component may be expected. As described elsewhereherein, techniques may be used in determining an HI in accordance with avector of CI values having expected CI values included in vector M_(i)for a given HI classification, i, having a covariance matrix Σ_(i). Onetechnique uses a three state Kalman filter for predicting or trendingfuture HI values.

The Kalman filter may be used for various reasons due to the particularfactors taken into account in the embodiment and uses described herein.It should be noted that other systems embodying concepts and techniquesdescribed herein may also take into account other noise factors. In oneembodiment, the Kalman filter may be preferred in that it provides fortaking into account the noise of a particular arrangement of components.There may be noise corruption, such as indicated, for example, by thecovariance matrices described and used herein. It may be desirous tofilter out such known noise, such as using the Kalman filter, providingfor smoothing of data values.

The Kalman filter provides a way to take into account other aprioriknowledge of the system described herein. In particular, the health of acomponent, for example, may not change quickly with time. The differencebetween the health of a component at a time t, and time t+delta may notbe large. This technique may also be used in connection with determiningfuture HIs of a particular part, for example, where the part is old. Apart may have reached a particular state of relatively bad health, butstill a working and functional part. The techniques described herein maybe used with an older part, for example, as well as a newer part.

In the arrangement with the Kalman filter, state reconstruction may beperformed using the Ricatti equation, as known to those of ordinaryskill in the art. The technique that is described herein uses athree-state Kalman filter of HI, and the first and second derivativesthereof with respect to changes in time, denoted, respectively, dt² anddt³. The Ricatti equation in this instance uses a [1×3] vector of timevalues rather than a single value, for example, as may be used inconnection with a single state Kalman filter.

What will now be described are equations and values that may be used indetermining a future value of a particular HI. Let: $H = \begin{bmatrix}1 & 0 & 0\end{bmatrix}$ $\Phi = \begin{bmatrix}1 & {dt} & \frac{{dt}^{2}}{2} \\0 & 1 & {dt} \\0 & 0 & 1\end{bmatrix}$ $Q = {2\sigma^{2}{\overset{\_}{t}\begin{bmatrix}\frac{{dt}^{3}}{2} & \frac{{dt}^{2}}{2} & \overset{\_}{t} \\\frac{{dt}^{2}}{2} & {dt} & 1 \\\overset{\_}{t} & 1 & \frac{1}{\overset{\_}{t}}\end{bmatrix}}}$ $X = \begin{bmatrix}{HI\_ est} \\{\overset{.}{H}I} \\{\overset{¨}{H}I}\end{bmatrix}$

in which:

σ is the power spectral density of the system,

R is the measurement error,

P is the covariance,

Q is the plant noise,

H is the measurement matrix,

K is the Kalman gain and

Φ is the state transition matrix.

H may be characterized as the Jacobian matrix. Since the value of asingle HI is desired, only the first entry in the H vector is 1 withremaining zeroes. There are n entries in the n×1 vector H for the nstate Kalman filter. Similarly, the X vector above is column vector of 3HI entries in accordance with the three-state Kalman filter. The endvalue being determined is the vector X, in this instance whichrepresents a series of HI values, for which the first entry, HI_est inthe vector X is the one of interest as a projected HI value beingdetermined Within the vector X, H{dot over (I)} represents the firstderivative of HI_est and HÏ represents the second derivative of HI_est.{overscore (t)} represents the average amount of time betweenmeasurements or updating of the HI value. In other words, if dtrepresents a measurement or delta value in units of time between HIdeterminations, and this is performed for several instances, {overscore(t)} represents the average of the delta values representing timechanges.

What will now be presented are equations representing the relationshipsbetween the above quantities as may be used in determining a value ofX(1) for predicting or estimating an HI value at a future point in timegiven a current HI value.

X _(t|t−1) =ΦX _(t−1|t−1)  (Equation T1)

P _(t|t−1) =ΦP _(t−1|t−1)Φ^(T) +Q  (Equation T2)

K=P _(t|t−1) H ^(T)(HP _(t|t−1) H ^(T) +R)  (Equation T3)

P _(t|t)=(I−KH)P _(t|t−1)  (Equation T4)

X _(t|t) =X _(t|t−1) +K(HI−HX _(t|t−1))  (Equation T5)

Note that the subscript notation above, for example, such as “t|t−1”refers to determining a value of at a time t conditioned on themeasurement at a time of “t−1”. Similarly, “t|t” refers to, for example,determining an estimate at a time “t” conditioned on a measurement oftime “t”.

The current HI determined, for example, using other techniques describedherein, may be input into Equation T5 to obtain a projected value forHI_est, the best estimate of the current HI. To project the expect HI“n” units of time into the future, input the number of units of time“dt” into Φ (as described above), and use the state update equation(Equation T1) where now Equation T1 becomes: X_(t+dt|t) =ΦX _(t|t). Thisallows the best prediction of HI_est any number of units of time intothe future where HI_est is desired. It should be noted that as set forthabove, the linear matrix operation such as Φ X is equivalent to anintegration from t to dt of the state of X, where X represents thevector of HI values set forth above.

Different values may be selected for initial conditions in accordancewith each embodiment. For example, an initial value for P representingthe covariance may be (1/mean time value between failures). Anembodiment may use any one of a variety of different techniques toselect an initial value for P. Additionally, since P converges rapidlyto an appropriate value and the time between data acquisitions is smallin comparison to the mean failure time, selecting a particularly goodinitial value for P may not be as important as other initial conditions.A value for σ may be selecting in accordance with apriori information,such as manufacturer's data on the mean time between component parts'expected failure time. For example, for a transmission, the mean failuretime may be approximately 20,000 hours. The spectral density may be setto (1/20,000)². It should be noted that the failure rates may begenerally characterized as an exponential type of distribution. The meantime between expected failures is a rate, and the variance is that rateto the second power. R may also be determined using apriori information,such as manufacturer's data, for example, an estimated HI variance ofmanufacturer's data of a healthy component. Q may be characterized asthe mean time between failures and dt (delta change in time betweenreadings). As the value of dt increases, Q increases by the third power.

Input data used in the foregoing trending equations may be retrievedfrom collected data, for example, as may be stored in the system of FIG.1.

In determining HIs, for example, as in connection with the system ofFIG. 1 for particular components, HIs may be derived using one or moreCIs. In calculating CIs, data acquisitions may occur by recordingobserved data values using sensors that monitor different components.There may be a need for estimating data used in connection with CIcalculations, for example, in instances in which there may be too littleor no observed empirical. For example, in connection with a power train,there may be a need to obtain estimated data, for example, for eachbearing, shaft and gear within the power train to calculate CIs.However, insufficient empirical data may exist in connection with gearor bearing related measurements, such as, for example, those inconnection with a gear or bearing related measurements, such as, forexample, those in connection with a gear or bearing fault due to therare occurrence of such events. In such instances, mean and thresholdvalues may be derived using other techniques.

A CI may indicate a level of transmission error, for example, in whichtransmission error is a measure of the change in gear rigidity andspacing. Modeling transmission error may allow one to gauge CIsensitivity and derive threshold and mean values indicative ofgear/bearing failure. This transmission error modeling may be referredto as dynamic analysis. What will now be described is a technique thatmay be used to model a gears to obtain such estimated values. Bymodeling each gear pair as a damped spring model with the contact linebetween the gears, transmission error may be estimated. It should benoted that this model uses two degrees of freedom or movement. Othersystems may use other models which may be more complex having moredegrees of freedom. However, for the purposes of estimating values, thismodel has proven accurate in obtaining estimates. Other embodiments mayuse other models in estimating values for use in a system such as thatof FIG. 1.

Referring now to FIG. 30, shown is an example of an illustration of apair of gears for which a model will now be described. A force P at thecontact gives linear and torsional response to each of the 2 gears for atotal of four responses as indicated in FIG. 30. The relative movement dat P is the sum of the 4 responses together with the contact deflectiondue to the contact stiffness s_(c) and the damping coefficient b_(c).This may be represented as: $\begin{matrix}{d = {P\left( {\frac{1}{{sp} + {{j\omega}\quad {bp}} - {m\quad p\quad \omega^{2}}} + \frac{{rp}^{2}}{{kp} + {{j\omega}\quad {qp}} - {{Ip}\quad \omega^{2}}} + \frac{1}{{sw} + {{j\omega}\quad {bw}} - {{mw}\quad \omega^{2}}} + \frac{{rw}^{2}}{{kw} + {j\quad \omega \quad p\quad w} - {{Iw}\quad \omega^{2}}} + \frac{1}{{sc} + {j\quad \omega \quad {bc}}}} \right)}} & \text{EQUATION~~G1}\end{matrix}$

in which:

sp is the linear stiffness of the pinion;

j is the square root of −1;

ω is the angular rate that may be obtained from the configuration file(e.g., shaft rpm*60*2π to obtain radians per second for the piniondriving the wheel);

bp is the linear damping coefficient of the pinion;

mp is the mass of the pinion;

rp is the radius of the pinion;

kp is the angular effective stiffness of the pinion;

qp is the angular damping coefficient of the pinion;

Ip is the angular effective mass of the pinion;

sw is the linear stiffness of the wheel;

bw is the linear damping coefficient of the wheel;

mw is the mass of the wheel;

rp is the radius of the pinion;

kw is the angular effective stiffness of the wheel;

qw is the angular damping coefficient of the wheel;

Iw is the angular effective mass of the wheel;

sc is the linear stiffness of the contact patch where the two gears comeinto contact;

bc is the linear damping coefficient of the contact patch;

It should be noted that values for the above-referenced variables on therights hand side of EQUATION G1 above, except for P (described below),may be obtained using manufacturer's specifications for a particulararrangement used in an embodiment. An embodiment may include quantitiesfor the above-referenced variables in units, for example, such asstiffness in units of force/distance (e.g., newtons/meter), mass in kgunits, and the like.

The relative movement, d, is the T.E., so from d, the above-referencedequation can be solved for P, the tooth force. Deflection is the force(input torque divided by the pinion base radius)*the elastic deflectionof the shafts, which may be used in estimating P represented as:

P=(1/kp*rp)+(1/sp)+(1/sw)+(1/kwrw)  EQUATION G2

where the variables are as described above in connection with EQUATIONG1. Using the above estimate for P with EQUATION G1, the displacement,such as a vibration transmitted through the bearing housing andtransmission case (which acts an additional transfer function), may bedetermined.

Referring again to FIG. 30, shown is an example of an illustration ofthe gear model and the different variables used in connection withEQUATION G1 and G2. Lp may represent the longitudinal stiffness of thepinion and Lw may represent the longitudinal stiffness of the wheel. Itshould be noted that these elements may not be included in an embodimentusing the two degrees of freedom model.

Bearings may also be modeled to obtain estimates of fault conditions ininstances where there is little or no empirical data available. Withbearings, a periodic impulse is of interest. The impulse is the resultof a bearing rolling over a pit or spall on the inner or outer bearingrace. The intensity of the impulse on the bearing surface is a functionof the angle relative to the fault, which may be represented as, forexample, described in the Stribeck equation in a book by T. A. Harris,1966, Rolling Bearing Analysis. New York: John Wiley p 148 as:

q(θ)=q ₀[1−({fraction (1/2ε)})(1−cos θ)]^(n)  EQUATION B1

where n=3/2 for ball bearings and 10/9 for rolling elements bearing,ε<0.5, and θ is less than π/2 in accordance with values specified inthis particular text for the different bearings used in theabove-referenced Stribeck equation represented as in EQUATION B1.

An impulse in a solid surface has an exponential decay constant, whichmay be taken into account, along with a periodic system due to rotationof the shaft. The bearing model may then be represented as a quantity,“s”, which is the multiplication of the impulse, “imp” below, theimpulse intensity, “q(θ)” as may be determined above, the period shaftrotation, which is “cos(θ)” below, all convoluted by the exponentialdecay of the material and represented as:

s=[imp×q(θ)×cos θ]exp(T/t)  EQUATION B2

where T is the exponential decay and t is the time. It should be notedthat “T” varies with the material of the solid surface. “exp(T/t)” maybe obtained, for example, using a modal hammer, to generate the decayresponse experimentally. An embodiment may also obtain this value usingother information as may be supplied in accordance with manufacturer'sinformation. The value of “t” may be a vector of times starting with thefirst time sample and extending to the end of the simulation. T isgenerally small, so the expression “exp(T/t)” approaches zero rapidlyeven using a high sampling rate.

“imp” is the impulse train that may be represented as the shaftrate*bearing frequncy ratio*sampling rate for the simulation period.

“s” is the simulated signal that may be used in determining a spectrum,“S”, where “S=fft(s)”, the Fourier transform of s into the frequencydomain from the time domain. As described in more detail in followingparagraphs, in determining a CI in connection with the bearing modelsignal “s” having spectrum “S”, for example, the Power Spectral Densityof S at a bearing passing frequency may be used as a CI. Additionally,for example, other CI values may be obtained, such as in connection withthe CI algorithm comparing the spectrum “S” to those associated withtransmission error in connection with a normal distribution using thePDF/CDF CI algorithms that may be generally described as hypothesistesting techniques providing a measure of difference with regard whetherthe spectrum is normally distributed.

It should be noted that, as described elsewhere herein in connectionwith gear models, values may be used in the foregoing equations inconnection with simulating various fault conditions and severity levels.The particular values may be determined in accordance with what smallamount of observed data or manufacturer's data may be available. Forexample, in accordance with observed values, an impulse value of 0.02for the impulse, “imp”, may correspond to a fairly severe faultcondition. Values ranging from 0.001 to 0.03, for example, may be usedto delimit the range of “imp” values used in simulations.

Following is an example of estimated data using the foregoing equationsfor a bearing having the following configurations:

Rpm=287.1

Roller diameter=0.25

Pitch diameter=1.4171

Contact angle=0

Number of elements=10

Inner race fault

Referring now to FIG. 31, shown is an example of a graphicalrepresentation of the signal for the foregoing configuration when thereis some type of bearing fault as estimated using the foregoing equationsEQUATION B1 and B2. FIG. 32 represents the estimated spectrum “S” as maybe determined using EQUATION B2 above.

It should be noted that for bearings, there may be three types offaults, for example, estimated using the foregoing equations. There maybe an inner race fault, an outer race fault or a roller element fault.Localized bearing faults induce an excitation which can be modeled as animpulse train, expressed as imp in the above equation. This impulse“imp” corresponds to the passing of the rolling elements of the fault.Assuming a constant inner ring rotation speed, the impulse train isperiodic and the periodicity depends on the fault location.

For outer race faults, the bearing frequency ratio, f_(d), or may berepresented as: $\begin{matrix}{f_{d,{or}} = {\frac{N}{2}\left( {1 - {\frac{d_{b}}{d_{m}}{\cos (\alpha)}}} \right)\left( {f_{ir} - f_{or}} \right)}} & \text{EQUATION~~B3}\end{matrix}$

where:

“d_(b)” represents the roller diameter,

“d_(m)” represents the pitch diameter,

“α”=2π*frequency, f_(d);

“f_(ir)” is the rotation frequency of the inner race (e.g. shaft rate),and

“f_(or)” the rotation frequency of the outer race (if fixed=0).

For inner race faults, the bearing frequency ratio, f_(d, ir) may berepresented as: $\begin{matrix}{f_{d,{ir}} = {\frac{N}{2}\left( {1 + {\frac{d_{b}}{d_{m}}{\cos (\alpha)}}} \right)\left( {f_{ir} - f_{or}} \right)}} & \text{EQUATION~~B4}\end{matrix}$

Replacing α with 2πf_(d), the time response is f(t). This substitutionmay be performed as the initial value of α is based on an angle and nota function of time. In a simulation, there is a time dependent responseas expressed using f(t).

The radial load applied to the bearing is not constant and results in aload distribution, which is a function of angular position. If thedefect is on the outer race, the amplitude of the impulse is constantbecause the fault location is not time varying. For an inner race fault,the amplitude with respect to angular position. The function is:$\begin{matrix}{{q(\theta)} = {{\begin{Bmatrix}\left( \left. {q_{o}\left( {1 - {\frac{1}{2ɛ}\left( {1 - {\cos \quad \theta}} \right)}} \right.} \right) \right)^{n} \\{0\quad \text{elsewhere}}\end{Bmatrix}\quad \text{for}\quad {\theta }} \leq \theta_{\max}}} & \text{EQUATION~~B5}\end{matrix}$

q(t)=q(2πf(θ)). This quantity q(t) is amplitude at a particular time, orq(theta) representing the amplitude at a particular angle. Amplitudemodulation takes into account the distance from the fault to the sensor.For outer race fault, the quantity cos (θ) is constant (1), for innerrace fault, it is the cosine function, noted as “cos (θ)” in the aboveequation.

For a linear system, the vibrations at a given frequency may bespecified by the amplitude and phase of the response and the timeconstant of the exponential decay. As the angle, θ above, changes, theimpulse response, h(t), and the transfer function H(f) also change dueto the changing transmission path and angle of the applied impulse. Itis assumed that the exponential decays is independent of the angle θ, sothat the response measured at a transducer due to an impluse applied tothe bearing at the location θ is characterized by an amplitude which isa function of θ.

The impulse response function h(t) and the transfer function H(f) may bereplaced by a function a(θ) giving the amplitude and sign of thetransfer function H(f)

at each angle theta and by the exponential decay of a unit impulse,(e(t)). For an inner race defect, rotating at the shaft frequency fs,the instantanous amplitude of the transfer function between the defectand the transducer as a function of time, a(t) may be obtained bysubstituting 2π*fs*t for theta. Note that a(t) is periodic. At θ=0relative to the defect and transducer, a(t) has its maximum value. Atθ=π, a(t) should be a minimum because the distance form the defect tothe transducer is a minimum. Additionally, the sign is negative becausethe impulse is in the opposite direction. Because of these properties,the cos(t) may be used for the function a(t).

The impulse train is exponentially decaying. The decay of a unit impulsecan be defined by:

e(t)=exp(−t/T _(e))  EQUATION B6

for t>0, where T_(e) is the time constant of decay.

The bearing fault model is then:

v(t)=[imp(t)q(t)a(t)]*e(t)  EQUATION B7

where:

imp(t), which is the impulse over a time t,=2π*shaft rate*time*bearingfrequency ratio, as may be determined using EQUATIONS B3 and B4 above;

a(t) is the cos(θ) for an inner race, which is 1 for an outer race,where cos(θ)=0, where θ is time varying;

and q(t) and e(t) are as described above.

An embodiment may include a signal associated at the sensor for gear andbearing noise combined from the bearing and the gear model may berepresented as:

s(t)=[d(t)f(t)q(t)a(t)]*e(t)*h(t)  EQUATION B8

where:

h(t) is the frequency response of the gear case, as may be determined,for example, using an estimate produced with linear predictive coding(LPC) techniques or with a modal hammer analysis;

d(t) is the signal associated with gear/shaft T.E. as may be determinedusing the gear model EQUATION G1;

and other variables are as described elsewhere herein.

The frequency spectrum of signals representing a combined bearing andgear model from EQUATION B8 may be represented as:

 S(f)=[D(f)*F(f)*Q(f)*A(f)]E(f)H(f)  EQUATION B9

As described elsewhere herein, healthy data, such as may be obtainedusing manufacturer's information, may be used in determining differentvalues, such as those in connection with stiffnesses for gearsimulation, amplitude and exponential decay for bearing faults. In termsof generating fault data, since these systems are linear, the followingmay be defined:

For gear faults indicative of a crack, a reduction in the stiffness fora tooth (e.g. 50 and 20 percent of normal) may be used in estimatingmedian and high fault values. Additionally, these values may be varied,for example, using the Monte Carlo simulation to quantify variance.

For shaft misalignment, shaft alignment within the model may be variedto estimate mean fault values

For gear spalling faults, the “size” of an impulse may be determinedthrough trial and error, and by comparing simulation values with anylimited observed fault data previously collected.

For bearing fault models, which are spalling faults, the size of animpulse, indicative of a fault, with known bearing faults, may bedetermined similarly as with gear spalling faults

Sensitivity analysis may be performed, for example using range ofdifferent input values for the different parameters, to provide forincreasing the effectiveness of fault detection techniques, for example,as described and used herein. For example, an embodiment may be betterable to simulate a family of bearing faults to tailor a particular CIalgorithm to be sensitive to that particular fault.

Using the foregoing, the modulated transmission error of a gear mesh,for example, which is a signal may be simulated or estimated. Thissignal may subsequently be processed using any one or more of a varietyof CI algorithms such that estimates for the mean and threshold valuescan then be derived for fault conditions. (It is assumed that thestiffness and torque are known apriori). Parameter values used in theabove equations corresponding to a healthy gear, for example, as may bespecified using manufacturer's data, may be modified to estimateparameter values in connection with different types of faults beingsimulated. By modifying these parameter values, different output valuesmay be determined corresponding to different fault conditions.

For example, known values for stiffness, masses, and the like used inEQUATION G1 may be varied. A cracked gear tooth may be simulated bymaking the stiffness time varying. The contact pitch may be varied withtime in simulating a shaft alignment fault. A modulated input pulse on dmay be used in simulating a spall on a gear tooth. Different parametervalues may be used in connection with specifying different degrees offault severity, such as alarm levels and warning levels. A particularparameter value, such as a tooth stiffness of 70% of the normalmanufacturer's specified stiffness, may be used in simulating warninglevels. A value of 20% of the normal manufacturer's specified stiffnessmay be used in simulating alarm levels. The particular values may bedetermined in accordance with comparing calculated values with thecharacteristics of real CI data on any few real faults collected.

While the invention has been disclosed in connection with the preferredembodiments shown and described in detail, various modifications andimprovements thereon will become readily apparent to those skilled inthe art. Accordingly, the spirit and scope of the present invention isto be limited only by the following claims.

What is claimed is:
 1. A method for normalizing at least one set of oneor more observed condition indicators comprising: determining aplurality of condition indicators and at least one relevant factorassociated with the condition indicators in accordance with previousdata acquisitions; determining a mean of said plurality of conditionindicators and at least one model coefficient corresponding to said atleast one relvant factor; and adjusting at least one new set of one ormore observed condition indicators in accordance with said at least onemodel coefficient and said at least one relevant factor to produce anormalized set of one or more indicators.
 2. The method of claim 1,wherein said at least one relevant factor includes at least one oftorque, airspeed, and rotor speed.
 3. The method of claim 2, furthercomprising: solving for B using the equation form y=B x for thefollowing: $y = {{\begin{bmatrix}{CI}_{1} \\{CI}_{\ldots} \\{CI}_{n}\end{bmatrix}\quad B} = {{\begin{bmatrix}B_{0} & {B_{1}\quad \ldots \quad B_{n}}\end{bmatrix}\quad \text{and}\quad x} = {\begin{bmatrix}1 & t_{1} & N_{R_{1}} & {Airspeed1} \\1 & {t\quad \ldots} & N_{R\quad \ldots} & {{Airspeed}\quad \ldots} \\1 & t_{n} & {N_{R}n} & {Airspeedn}\end{bmatrix}}}}$

wherein the vector y includes condition indicator values and the matrixx includes said at least one relevant factor corresponding to saidcondition indicator values included in y, and B includes coefficients.4. The method of claim 3, wherein a normalized condition indicator (CI)is determined at a later point in time by determining:CInormalized=CIobs−(B*x) where CIobs is an actual observed valuecorresponding to a normalized CI value CInormalized.
 5. The method ofclaim 4, further comprising: storing the values for B and x in adatabase; and retrieving said values for B and x in connection with saiddetermining of the CInormalized.
 6. The method of claim 3, wherein B isdetermined using x and y as follows: B=(x ^(T) x)⁻¹ x ^(T) y and x^(T)is a transpose of x.
 7. The method of claim 3, wherein B is determinedby minimizing an error between said plurality of condition indicatorsand a set of estimated condition indicator values or model conditionindicator values.
 8. The method of claim 7, wherein said error is a sumof square error between said plurality of condition indicators and saidset of estimated condition indicator values or model condition indicatorvalues.
 9. The method of claim 3, wherein B₀ is said mean.
 10. Acomputer program product for normalizing at least one set of one or moreobserved condition indicators comprising machine executable codes for:determining a plurality of condition indicators and at least onerelevant factor associated with the condition indicators in accordancewith previous data acquisitions; determining a mean of said plurality ofcondition indicators and at least one model coefficient corresponding tosaid at least one relevant factor; and adjusting at least one new set ofone or more observed condition indicators in accordance with said atleast one model coefficient and said at least one relevant factor toproduce a normalized set of one or more condition indicators.
 11. Thecomputer program product of claim 10, wherein said at least one relevantfactor includes at least one of torque, airspeed, and rotor speed. 12.The computer program product of claim 11 further comprising machineexecutable codes for: solving for B using the equation form y=B x forthe following: $y = {{\begin{bmatrix}{CI}_{1} \\{CI}_{\ldots} \\{CI}_{n}\end{bmatrix}\quad B} = {{\begin{bmatrix}B_{0} & {B_{1}\quad \ldots \quad B_{n}}\end{bmatrix}\quad \text{and}\quad x} = {\begin{bmatrix}1 & t_{1} & N_{R_{1}} & {Airspeed1} \\1 & {t\quad \ldots} & N_{R\quad \ldots} & {{Airspeed}\quad \ldots} \\1 & t_{n} & {N_{R}n} & {Airspeedn}\end{bmatrix}}}}$

wherein the vector y includes condition indicator values and the matrixx includes said at least one relevant factor corresponding to saidcondition indicator values included in y, and B includes coefficients.13. The computer program product of claim 12, wherein a normalizedcondition indicator (CI) is determined at a later point in time, and thecomputer program product further comprising machine executable codes fordetermining: CInormalized=CIobs−(B*x) where CIobs is an actual observedvalue corresponding to a normalized CI value CInormalized.
 14. Thecomputer program product of claim 13, further comprising machineexecutable codes for: storing the values for B and x in a database; andretrieving said values for B and x in connection with said determiningof the CInormalized.
 15. The computer program product of claim 12,further comprising machine executable codes for determining B using xand y as follows: B=(x ^(T) x)⁻¹ x ^(T) y and x^(T) is a transpose of x.16. The computer program product of claim 12, further comprising machineexecutable codes that determine B by minimizing an error between saidplurality of condition indicators and a set of estimated conditionindicator values or model condition indicator values.
 17. The computerprogram product of claim 16, wherein said error is a sum of square errorbetween said plurality of condition indicators and said set of modelcondition indicator values or estimated condition indicator values. 18.The computer program product of claim 12, wherein B₀ is said mean.